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FOREWORD 


This  report  describes  results  of  research  performed  under  Contract  No.  F33615-86-C-3612 
entitled  "Nonlinear  Flying  Qualities,"  prepared  for  the  Flying  Qualities  Group  of  the  Control 
Dynamics  Branch  (FIGCB)  of  the  Air  Force  Wright  Aeronautical  Laboratories  (AFWAL).  The 
general  objective  of  the  program  was  to  develop  and  evaluate  analytical  methods  which  can  be 
used  to  calculate  Nonlinear  Flying  Quality  Parameters  (NLFQFs).  The  specific  objective  was  to 
provide  interactive  computer-aided  analysis  tools  (based  on  the  nonlinear  inversion  concept)  for 
evaluating  nonlinear  flying  qualities  of  current  airplanes  and  guiding  the  development  of  future 
airplanes  with  improved  flying  qualities. 

The  research  was  performed  at  the  Honeywell  Systems  and  Research  Center  between  30  June 
1986  and  30  June  1987.  At  Honeywell,  the  research  was  conducted  by  Drs.  C.A.  Harvey  and 
B.G.  Morton,  Principal  Investigators;  Mr.  M.R.  Elgersma  and  Ms.  G.  Hines,  Associate 
Investigators;  and  Dr.  M.F.  Barrett,  Program  Manager.  Dr.  G.R.  Sell,  a  professor  in  the  School 
of  Mathematics  of  the  University  of  Minnesota,  served  as  a  consultant.  The  AFWAL/FIGCB 
technical  administration  was  performed  by  Mr.  C.F.  Suchomel,  Project  Engineer. 
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SECTION  1:  INTRODUCTION 


The  usefulness  of  current  flying  qualities  parameters  is  limited  by  their  almost  complete 
dependence  on  linear  analysis.  This  dependence  is  based  on  mathematical  tractability  and 
analysis  of  a  class  of  small  amplitude  maneuvers  about  trim  conditions.  The  limitations  of 
existing  flying  quality  parameters  leave  expensive  and  extensive  flight  testing  as  the  only 
accurate  means  of  assessing  flying  quality  during  highly  dynamic  maneuvers.  Inadequacies 
found  during  flight  test  can  lead  to  extremely  costly  modifications  and  redesigns.  This 
motivated  our  program  aimed  at  the  development  of  flying  quality  parameters  which  take  into 
account  the  nonlinearities  encountered  during  current  and  future  combat  maneuvers. 

We  took  a  novel  approach  to  provide  interactive  and  computer-aided  analysis  tools  which 
can  be  used  in  developing  and  evaluating  nonlinear  flying  quality  parameters  (NLFQPs).  This 
approach  is  based  on  a  technology  called  dynamic  inversion  which  we  used  to  generate 
maneuvers  and  analyze  flying  qualities.  To  enhance  our  analytical  developments  and  to 
expose  the  issues  concerning  NLFQPs  to  the  academic  community.  Professor  George  R.  Sell 
served  as  a  consultant. 

Our  technical  program  consisted  of:  a  review  of  existing  techniques,  problem  formulation, 
technical  development,  validation  and  illustration,  hi  the  problem  formulation,  the  analytical 
structure  of  models  was  defined  along  with  preliminary  definitions  of  sets  of  maneuvers  and 
potential  NLFQPs  to  be  examined.  The  technical  approach  consisted  of  analytic,  algorithmic, 
and  software  development.  Maneuvers  were  flown  in  simulation  during  the  validation  and 
illustration  to  demonstrate  the  utility  of  the  technical  development. 

The  concept  of  dynamic  inversion  is  quite  simple.  Suppose  the  aircraft  model  is  defined 
by: 


x  =  f(x)  +  6 


where  f  represents  the  uncontrolled  aircraft  dynamics  and  5  is  the  commanded  actuator  signal. 
For  a  desired  aircraft  response: 


x  =  L(x)  +  (pilot  command) 


dynamic  inversion  is  a  control  structure  which  forces  the  aircraft  to  respond  as  desired. 
Dynamic  inversion  is  done  as  follows: 

1)  Measure  the  state  x 

2)  Compute  f(x)  -  L(x) 

3)  Generate  the  actuator  command  signal:  5  =  L(x)  -  f(x)  +  (pilot  command) 


For  details,  see  section  5.1  . 
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SECTION  2:  OVERVIEW  AND  SUMMARY  OF  REPORT 


This  report  documents  the  results  of  our  efforts  to  develop  tools  that  can  be  used  for  the 
computation  of  nonlinear  flying  quality  parameters.  We  have  developed  many  new  ideas  for 
approaching  this  problem.  In  the  course  of  our  activities,  we  have  found  a  variety  of  candidate 
nonlinear  flying  quality  parameters  and  candidate  specifications  for  them.  These  parameters 
are  genuinely  different  from  expressions  derived  from  linearized  models:  we  work  with  the 
nonlinear  aerodynamic  functions  themselves  and  not  their  derivatives.  We  believe  the  candi¬ 
date  specifications  we  have  outlined  in  this  report  could  be  applied  to  current  and  future  Air 
Force  vehicles  to  improve  their  flying  qualities.  The  parameters  we  have  defined  here  can  be 
computed  directly  from  preliminary  nonlinear  aircraft  models. 

First  we  give  a  brief  outline  of  the  sections  to  come,  and  then  follow  with  a  more 
detailed  outline.  Section  3  shows  the  nonlinear  aircraft  models  we  used.  Sections  4  and  S 
present  the  theoretical  techniques  we  used  while  working  with  the  nonlinear  models.  Most  of 
these  techniques  were  developed  during  the  course  of  the  program  to  help  us  identify  and 
compute  the  parameters  discussed  in  section  6.  Anyone  interested  in  getting  to  the  flying  qual¬ 
ities  right  away  can  go  straight  to  section  6  and  look  through  sections  3,  4,  and  3  as  needed. 
Section  7  contains  the  simulation  results.  Section  7  is  followed  by  the  reference  list,  and  by 
the  appendices  treating  trajectories  and  dynamical  properties  of  maneuvers. 

Section  3  is  a  complete  description  of  these  models  in  the  form  used.  These  are  the  stan¬ 
dard  equations  of  motion  in  a  form  which  is  not  completely  general,  but  general  enough  to 
capture  most  of  the  important  features  of  the  models  commonly  in  use.  Our  techniques  will 
apply  to  fully  general  models,  but  only  at  the  expense  of  increased  analytical  complexity.  A 
description  of  the  application  of  dynamic  inversion  to  the  general  models  can  be  found  in  sec¬ 
tion  3.4  of  the  technical  proposal  for  this  contract  (in  response  to  PRDA  86-1  PMRN). 

The  one  special  assumption  made  is  that  the  effect  of  the  m  control  inputs  to  the  aircraft 
can  be  represented  by  a  6-by-m  matrix  function  of  the  state  that  multiplies  an  m  dimensional 
vector  function  of  the  state  and  the  control  input  signals.  This  is  the  situation  that  arises  when, 
for  example,  the  controls  are  described  by  aerodynamic  derivatives  such  as  C,^  that  are  func¬ 
tions  of  the  state. 
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Section  4  discusses  a  novel  approach  to  the  computation  of  trim  conditions  for  the 
models  discussed  in  section  3.  Here  we  exploit  the  special  form  of  the  equations  to  reduce  the 
computations  to  an  algorithm  that  is  simple  enough  to  allow  analytic  computations.  Sections 
4.1  and  4.2  explain  the  basic  ideas  behind  this  approach  and  present  the  class  of  models  to 
which  they  can  be  applied.  Section  ^.3  makes  the  algorithm  explicit  for  a  family  of  systems 
that  includes  those  presented  in  section  3.  Section  4.4  explains  the  two  different  ways  that  we 
have  implemented  these  ideas  for  the  solution  of  aircraft  equilibria.  Section  4.5  gives  the 
details  of  the  computations  for  the  aircraft  models,  using  both  approaches  discussed  in  section 
4.4.  Included  in  section  4.5  is  a  special  treatment  for  aircraft  exhibiting  leftAight  symmetry 
and  other  special  features  that  are  often  assumed  and  which  simplify  the  results. 

Section  5  is  our  analysis  of  nonlinear  aircraft  dynamics  using  partial  dynamic  inversion. 
Section  5.1  is  a  short  summary  of  the  method  of  partial  dynamic  inversion  and  five  examples 
of  inversion  approaches  that  we  investigated  on  this  program.  Section  5.2  introduces  the 
theory  of  complementary  dynamics  in  a  general  setting,  then  section  5.3  shows  the  algorithm 
for  computing  the  complementary  dynamics  for  the  nonlinear  airplane  models.  The  dimension¬ 
less  coordinates  are  introduced  here  to  make  it  possible  to  fit  the  equations  into  a  reasonable 
amount  of  space  and  to  reduce  the  number  of  independent  parameters  appearing  in  the  equa¬ 
tions  of  motion.  Then,  for  the  special  cases  mentioned  before,  the  special  form  of  the  com¬ 
plementary  dynamic  parameters  is  derived  explicitly.  It  is  shown  then  that  for  typical  models 
in  this  form,  the  a,P,0,<I>  complementary  dynamic  equations  are  stable  and  have  unique 
equilibria.  Even  if  the  uncontrolled  aircraft  is  unstable,  the  nonlinear  inverter  stabilizes  the  air¬ 
craft.  Section  5.4  then  shows  how  to  transform  the  general  a,P,@,d>  complementary  dynamic 
equations  into  canonical  forms,  so  that  Liapunov  functions  can  be  applied.  Section  5.5  then 
demonstrates  how  the  Liapunov  functions  are  constructed,  and  explains  the  significance  of  the 
result  to  aircraft  stability. 

The  theory  developed  in  sections  3,  4,  and  5  is  central  to  the  discussion  of  nonlinear 
flying  qualities  in  section  6,  but  it  is  not  a  prerequisite.  The  reader  is  warned  that  the  parame¬ 
ters  developed  in  section  6  are  represented  in  terms  of  the  dimensionless  coordinates 
developed  in  section  5.3,  and  that  conversions  to  account  for  units  in  the  final  answers  might 
be  necessary  before  meaningful  comparisons  can  be  made  between  different  aircraft. 

In  section  6.1  we  introduce  the  notion  of  commanded  dynamics,  then  present  a  list  of  the 
sections  from  the  MIL-F-8785C  document  that  we  believe  are  specifications  for  commanded 


dynamics.  Then  we  present  some  parameters  that  we  call  commanded  dynamic  parameters, 
derived  from  the  form  of  the  nonlinear  models.  These  include  the  basic  Q  command 
effectiveness  parameter,  the  dynamic  pitch-control  ratio,  the  minimum  lateral-directional  com¬ 
mand  effectiveness  parameter,  and  the  dynamic  lateral-directional  control  ratio  vector.  In  sec¬ 
tion  6.2  we  introduce  the  notion  of  complementary  dynamics,  then  present  a  list  of  the  sec¬ 
tions  from  the  MIL-F-8785C  document  that  we  believe  are  specifications  for  complementary 
dynamics.  Then  we  discuss  our  complementary  parameters  that  were  developed  in  complete 
detail  in  section  5.  Specific  comparisons  are  indicated  between  our  complementary  dynamic 
parameters  and  the  specifications  in  sections  3.2. 1.1  (longitudinal  static  stability),  3.4.1 
(dangerous  flight  conditions)  and  3.4.2  (flight  at  high  angle  of  attack)  of  MIL-F-8785C.  The 
discussions  at  the  end  of  cases  1  and  2  of  section  3.5  should  be  read  in  conjunction  with  sec¬ 
tion  6.2.  Section  6.3  is  a  short  but  very  interesting  discussion  of  lift-to-drag  ratios,  what  they 
look  like  for  the  F-4,  the  F-14,  and  the  F-15  aircraft  for  all  angles  of  attack,  and  why  that  has 
a  bearing  on  flying  quality.  Section  6.4  contains  some  criteria  for  coordinated  flight  at  high 
angle  of  attack  and  for  sustained  high-angle-of-attack  maneuvering.  In  section  6.5  we  present 
an  approach  towards  defining  dynamic  flying  quality  metrics  using  the  coordinated-flight 
U,P,Q,R  dynamic-inversion  controller  discussed  in  section  6.4  (A  prototype  of  this  controller 
is  demonstrated  in  the  simulations  of  sections  7.2  and  7.3).  In  section  6.6,  some  parameters 
for  measuring  the  influence  of  the  dynamic  aerodynamic  coefficients  on  flying  qualities  are 
given. 

Section  7  presents  some  maneuvers  generated  by  the  batch  simulation  using  nonlinear, 
partial  dynamic  inversion  controllers.  Section  7.1  contains  some  examples  of  a  roll  reversal, 
section  7.2  is  a  barrel  roll,  and  section  7.3  is  a  highly-dynamic  diving  turn.  We  list  the  main 
developments  that  arose  from  analysis  of  the  simulations  in  the  summary  section  7.4 

There  is  a  list  of  references,  and  then  the  two  appendices.  The  first  appendix  discusses 
the  trajectories  of  vehicles  at  equilibrium.  The  second  appendix  was  written  by  George  Sell. 
In  it,  he  proves  a  fundamental  lemma  for  the  general  theory  of  dynamic  inversion. 


SECTION  3  :  NONLINEAR  MODELS 


l 

i 


In  this  section  we  describe  the  aircraft  equations  of  motion.  For  more  details  on  the  nota¬ 
tion  and  derivation  of  the  aircraft  equations  of  motion,  see  [E3]. 

3.1  Notation 


Coordinates: 


U 

V 

W. 

rvi 


velocity  vector  of  die  c.g.  in  body-axis  coordinates 


P 

.a. 


M  = 

fPl 


=  velocity  vector  of  the  c.g.  in  wind-axis  coordinates 
a  (speed,  sideslip  angle,  angle  of  attack) 


Mach  number 


speed 

speed  of  sound 


BJ- 

'F.e/D 


angular  velocity  vector  in  body-axis  coordinates 


a  Euler  angles  for  heading,  elevation,  and  bank  angle 


(yaw,  pitch,  roll  sequence) 


T 

engine  thrust 

8. 

total  aileron  angle 

8« 

total  elevator  angle 

8r 

total  rudder  angle 

Tabular  aero  data  (from  wind  tunnel  testing): 


In  general  these  are  nonlinear  tabular  functions  of  many  variables,  e.g.  the  aero  force  in  the  x 
direction  is  given  by  Nx(M,V,a,p,d,p,P,Q,R,T,8t,8e,6r,  •  •  •  )  .  In  this  report,  we  have  only 
kept  the  M,V,a,P,P,Q,R,T,8lt8e,8r  dependence.  Furthermore,  we  have  expanded  the  functions 
in  a  Taylor  series  with  respect  to  P,Q,R,T,8t,8e,8r  and  kept  only  the  two  lowest  order  terms  in 
the  Taylor  series.  For  example, 

Nx(M,V,a,(5,P,Q,R,T,8t,8e,8r)  = 

V*pV2S  [Cx(M,a,P)  +  Cv(M,a,p)  P  +  C^fM.ct.p)  Q  +  QJM.a.p)  R  ]  + 

CXT(M,a,p)  T  +  VipV2S  [C^CM.a.P)  8,  +  C^M.a.p)  8e  +  C^M.a.P)  8r  ]  + 

higher  order  terms. 

When  linear  analysis  is  done,  the  above  types  of  expressions  are  expanded  further  in  a  Taylor 
series  with  respect  to  the  velocity  components  too  (e.g.  M,  a,  P  or  U,  V,  W  )  and  only  the 
two  lowest  order  terms  kept.  In  this  report,  we  will  be  working  directly  with  the  nonlinear 
functions  of  M,  a  ,  and  P  instead  of  expanding  them  in  a  Taylor  series  in  M,  a  ,  P  . 


Cx(M,a,P) 

Cy(M,a,P) 

C^M.a.P) 


=  static  nonlinear  aero  force  functions  in  body-axis  coordinates. 


C,(M.o,P) 

Cm(M,a,p) 

Cn(M,oc,p) 


=  static  nonlinear  aero  moment  functions  in  body  axes 


Note:  the  above  nonlinear  functions  of  M,  a  ,  and  P  are  the  coefficients  on  the  zeroeth  order 
terms  in  the  Taylor  expansion  with  respect  to  P,Q,R,T,8t,8e,8r 


Cx,(M,o,P) 

Cyp(M,a,P) 

C,„(M,a,P) 


CXo(M,a,P) 

Cy„(M,a,P) 

ClQ(M,a,P) 


C«ll(M,<x,p) 

Cyx(M.a,p) 

Qp(M,a,p) 


=  dynamic  nonlinear  aero  force  functions  in  body 


Cv(M,a,P)  C^M.a.P)  C|JI(M,cx,P) 
C^M.a.P)  Cm,(M,a,P) 
Cnr(M,a,p)  C^CM.a.p)  C^M^p) 
axes 


dynamic  nonlinear  aero  moment  functions  in  body 


CXt(M,<x,P)  C^(M,a,p)  C^(M,a.P)  C^fM.a.P) 

CyM.a.P)  Cyn(M,a,P)  Cy^(M,a,P)  Cy.CM.a.p)  = 

(^(MaP)  C74(M,a,p)  CZ4(M,a,P)  C^(M,a,p) 

•  « 

nonlinear  aero  coefficients  describing  forces  due  to  the  controls  in  body  axes 


CIt(M.o,P)  Cu  (M,a,p)  Cu>(M,a,P)  Cu(M,a,p) 
Cmr(M.a,P)  Cnu  (M,a,P)  CIBi<(M,a,p)  C„^(MaP)  = 
Cnr(M,a,P)  Cni(M,a,p)  C^M.a.p) 

k  « 

nonlinear  aero  coefficients  describing  moments  due  to  the  controls  in  body  axes 


Note:  the  above  nonlinear  functions  of  M,  a  ,  and  p  are  the  coefficients  on  the  first  order 
terms  in  the  Taylor  expansion  with  respect  to  P,Q,R,T,8t,8e,8r 


Physical  parameters: 


p  =  air  density 
g  =  gravity 

Physical  parameters  for  the  aircraft: 

m,c,b,S  a  mass,  mean  aerodynamic  chord,  wing  span,  wing  area 

^xx  ~^xy  — ^xz 

-Ixy  Iyy  -Iyj.  a  moment  of  inertia  matrix  in  body  axes 

— ^XZ  ~Iyz  Izz 
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3.2  Introduction  to  the  Equations  of  Motion 

Consider  a  6DOF  nonlinear  aircraft  model  having  four  control  inputs.  The  equations  of 
motion  are  often  given  in  a  mixed  system  of  coordinates.  There  are  eight  states  that  the 
forces  and  moments  depend  on.  The  rigid-body  mechanics  are  most  easily  expressed  in  the 
body-axis  velocity  (U,V,W)  and  angular  rate  (P,Q,R)  coordinates,  and  two  of  the  Euler  angles: 
bank  angle  and  elevation  (0,6)  .  The  forces  and  moments  do  not  depend  on  the  heading 
angle,  'F,  so  it  is  not  included  as  a  state.  The  result  of  not  including  ¥  as  a  state  is  that  the 
equilibria  will  generally  consist  of  vertical  helices  instead  of  just  straight  lines  in  inertial  space 
(see  appendix  A).  The  aerodynamics  are  more  easily  expressed  in  terms  of  wind-axis  velocity 
(V,0,a)  in  place  of  (U,V,W).  In  order  to  fit  the  aircraft  model  into  the  framework  presented 
earlier,  a  single  system  of  eight  states  will  be  used,  and  the  angular  rates  will  be  expressed  as 
derivatives  of  the  Euler  angles  (<b,0,'F)  in  place  of  (P,Q,R).  In  this  mixed  coordinate  system, 
the  equations  of  motion  of  an  aircraft  are  of  the  form  shown  in  Figure  3.1. 


Mixed  Coordinate 
Aircraft 

Integrator  Chains 


Figure  3.1  Equations  of  Motion  in  Mixed  Coordinate  System 
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In  the  wind-axis  and  Euler  angle  coordinates,  the  state  vector  is  x  =  (V,P,a,<I>,0,4/,<I>,0).  The 
form  of  the  equations  of  motion  in  this  case  are  shown  in  Figure  3.2 


Wind-Axis,  Euler-Angle 
Aircraft 

Integrator  Chains 


Figure  3.2  Equations  of  Motion  in  Wind-Axis  and  Euler  Angle  Coordinate  System 


3 3  Equations  of  Motion 


The  aircraft  equations  of  motion  [E3]  may  be  written  in  the  form 

0  =  f(x,x)  +  g(x)  h(x,u) 

where 

x  is  (V,|5,a,<i>,0,4',<t>,©)T  ,  u  is 


(3.3.1) 


T 

5. 

5e 

Sr 


f  is  6  x  1 


g  is  6  x  4  and  h  is  4  x  1 


All  the  tabular  aero  functions  which  are  denoted  by  capital  C’s  with  subscripts,  like  Cx  ,  , 

Q,  ,  etc.  axe  nonlinear  functions  of  (M,a,P)  .  M  is  a  function  of  V  and  air  temperature. 

The  top  three  rows  of  (3.3.1)  are  the  rigid-body  force  equations,  while  the  bottom  three  rows 

y  moment  equations.  The  quantity  f(x,x)  +  g(x)  h(x,u)  contains  the 

which  are  given  in  terms  of  the  state  x  by 


of  (3.3.1)  are  the  rigid-bod 

P 

U 

expressions 

Q 

and 

V 

IrJ 

LwJ 

P 

Q 

IKJ 


1  0  -sin(©) 

0  cos(<D)  cos(0)sin(<l>) 
0  — sin(<l>)  cos(9)cos(<l>) 


6 

e 

(3.3.2) 


»  • 

u 

cos(P)cos(a) 

V 

.w. 

=  V 

sin(p) 

.cos(P)sin(a) . 

(3.3.3) 


The  four  components  of  h(x,u)  are  the  forces  (divided  by  mg)  produced  by  the  controls 
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h(x,u)  = 


VfrpV2  S 


T 

mg 


mg 

KpV2  S 
mg 

V4pV2S 

mg 


sin(8«- 


sin(8e-<x) 

sin(8r-P) 


(3.3.4) 


This  form  for  h  is  an  example  of  the  way  that  h  can  depend  on  x  and  u  in  a  nonlinear  way. 
The  sin  functions  were  chosen  to  account  for  the  fact  that  surface  deflections  are  bounded 
even  if  the  surfaces  could  rotate  through  a  full  2n  radians.  The  shift  in  the  arguments  of  the 
sin  functions  account  for  the  change  in  airflow  at  the  surfaces  caused  by  the  aircraft’s  rolling 
motion  and  angle  of  attack.  The  h  vector  gets  premultiplied  by  the  6  by  4  g(x)  matrix,  so 
all  four  controls  can  influence  all  six  degrees  of  freedom. 


The  derivative  in  the  rotating  accelerating  frame  of  the  aircraft  (sometimes  referred  to  as  the 
covariant  derivative)  of  the  linear  momentum  is 


1 

0 

0 

A 

OR  -Q 

f"  « 

u 

■ 

0 

1 

0 

Q 

At 

-R  0  P 

m 

V 

k 

0 

0 

1. 

at 

© 

a- 

O' 

fc 

LwJ 

4 

T  -sin(0) 

-  mg  cos(0)sin(<&) 
lcos(0)cos(<I>) 


The  aerodynamic  forces  with  zero  control  input  (divided  by  mg)  are: 


l*pV2S 

mg 


C, 

Cy 

C, 


1 

2V 


Cxp 

^yr 

S  C*  C,j 


b  0  0 
0  c  0 
0  0b 


P 

IR 


So  the  top  three  rows  of  f(x,x)  are 


-1 

mg 


* 

• 

1 

0 

0 

H 

0  R  -Q 

■ 

U 

-sin(0) 

0 

1 

0 

u 

1 

-R  0  P 

m 

V 

-  mg 

cos(9)sin(d>) 

. 

k 

p 

0 

1. 

dt 

Q  -P  0  . 

4 

k 

w 

4 

,cos(9)cos(<t>) . 

4 

V4pV2S 

mg 


C,  ^  b  0  0  Tpl 

+  2V  ^  ^y»  ®  c  0  Q 

fcJ  [c^  C,,  C^,  Ip  o  »J  1*J 


(3.33) 


Note  that  no  p  dynamics  have  been  included.  However,  the  controller  (introduced  in  sec¬ 
tion  S.l)  can  be  assumed  to  be  using  a  measured  (so  varying)  value  of  p  in  these  equations. 
This  corresponds  to  treating  p  as  a  slowly  varying  parameter,  compared  to  the  rapidly  varying 
state.  A  consequence  is  that  the  equilibrium  helices  (see  Appendix  A)  will  vary  slowly  as 
altitude  is  lost. 

The  forces  (divided  by  mg)  produced  by  the  controls  give  the  top  three  rows  of  g(x)  multi¬ 
plied  by  h(x,u)  : 


Qlr 

<Vt  C*.  Cy..  Cy., 
CL 


mg 

sin(8t-  -S.) 
mg  ‘  2V  ’ 

sin(8e-a) 

mg 

My-*  sin(5r-{3) 
mg 


The  resulting  linear  momentum  equations  give  the  top  three  rows  of  0  =  f(x,x)  +  g(x)  h(x,u) 


1  0  0 


OIO7  -  -R  0 

p  0  ijdt  Iq  -p 


R  -Q 
0  P 
-P  0 


U  -sin(9) 

mV  -  mg  cos(9)sin(d>) 
.  LW  J  J  lcos(9)cos(<t>) 


'/toV2S 


lb  0  0 


Fy  +  ->v  rVr  Cy. 


L C*c  c**. 


0  c  0 
0  0  bl 


y* 

[i 


8 


I 

1 


V. 

» 

tfi 

§ 

% 

* 

9 
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FnDnimmcTjrwTi 


CxT  CU  CU 


Cyr  Cy^  Cyt'  Cy, 

^*r  Cj, 


mg 

*PV*  S  sin®.-  — ) 
mg  *  2V ' 

^eyis 

mg 

sm(S,-P) 

mg  _ 


(3.3.6) 


The  derivative  in  the  rotating  accelerating  frame  of  the  aircraft  (sometimes  referred  to  as  the 
covariant  derivative)  of  the  angular  momentum  is 


1  0  0 


0  1  0  j-  -R  0  P 

oo  ijdt  LQ  -p  0 


R  -Q 
0  P 


^xx  “^xy  ~^x2  fp 


-Ixy  Iyy  -Iyx  Q 


Ho  -Iyr  1= 


The  aerodynamic  moments  with  zero  control  inputs,  multiplied  by 


b  0  0 
0  c  0 
0  0b 


V4pV2S 

mg 


On  +  On,  Qtv, 


Qiq 


b  0  0 
0  c  0 
0  0b 


So  the  bottom  three  rows  of  f(x,x)  are 


b  0  0 
—  0  c  0 
mg  0  0  b 


1  0  0 


R  -Q 


0  1  0  j-  -R  0  P 

0  0  1 J dt  [q  -p  0 


[ 

*xy 

-Ixz 

XX 

!*y 

!yy 

I*z 

-V 

Izz 

(3.3.7) 
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■  .  '  ►  P  K.' 

-‘.y«yv 

■f. 


■  •  '  jk  ■ 

V  V  V  . 


.  •  >  •  „  »  . 


v  vv 


vV' 


‘Il  4..  I',  I  I  I 


V*pV2S 

mg 


C|1  ^  C|q  Ck 

On  +  2V  ^'m« 

■^*  ?*r  ^V}  ^ . 


b  0  0 
0  c  0 
0  0b 


The  moments  produced  by  the  controls,  multiplied  by 

b  0  0  -l 
-OcO 

mg  |p  0  b 

give  the  bottom  three  rows  of  g(x)  multiplied  by  h(x,u) : 


Ch  CU  CU.  CUr 
Q»r  ^>*1,  ^ 
Ql,  Q»4  Ql| 


mg 

mg  2V 

sin(8  -a, 
mg 

SN8.-P, 

mg 


So  the  angular  momentum  equations,  multiplied  by 

[b  0  0]"1 
—  0  c  0 
mg  0  0  b 


b  0  0 
0  c  0 
0  0b 


1  0  0 


R  -Q 


0  1  0  -R  0  P 

00  1 J dt  Lq  -p  0 


XX 

*yy 

_Iyz 

Kz 

-Iyr 

■/:pV2S 

mg 


+  — —  C  C  C 

2y  mP  ^rrvj  '“m« 


IS  ^n<3 


b  0  0 
0  c  0 
0  0b 


! 

f 

I 

* 

s 

$ 

1 


y 

fc 

p*f. 


% 


tx 


KW 


(3.3.8) 


Qi»t 

Q*r  Qi*  ^4 


*•  * 

sin(6e-a) 

mg 

sinft-u, 

mg 


r/i  -  fexpression(3.3.5)l 
'  ”  lpxpression(3.3.7)J 


(3.3.9) 


g(x)  = 


(3.3.10) 


Qt  ^*4,  ^*4,  ^Ur 

^Vr  Cy»t  ^y*.  ^y* 

C*T  ^z«t  *~*4.  ^i«t 
C.T  <\  cu>  cUr 
Qnr  Qrvii  Qm*[ 

Q»t  ^iii  ^n»  Qt i 


The  0  =  f  +  g  h  equations  represent  the  dynamics  while  the  kinematic  equations  for  the 
relationship  between  the  Euler  angle  rates  and  the  body-axis  rates  are  given  by  (3.3.2)  . 


SECTION  4  :  EQUILIBRIUM  MANIFOLD 


4.1  Introduction 

Systems  of  ordinary  differential  equations  (ODE)  with  no  control  inputs  typically  have  a 
set  of  isolated  equilibrium  points  (so  that  the  set  of  equilibrium  points  is  zero-dimensional)  . 
If  a  single  control  parameter  is  added,  then  for  each  fixed  value  of  the  control  parameter,  there 
will  be  a  corresponding  set  of  isolated  equilibrium  points.  The  family  of  equilibrium  points 
generated  in  this  way  forms  a  one-dimensional  set  Similarly,  systems  with  m  control  parame¬ 
ters  will  typically  have  m-dimensional  equilibrium  sets.  When  the  ODE  are  not  continuous 
with  respect  to  the  controls  or  other  parameters,  the  equilibrium  set  can  also  be  discontinuous. 
However,  the  equilibrium  equations  used  in  this  report  will  still  remain  valid,  they  just  change 
value  suddenly  when  the  discontinuity  is  reached.  Under  fairly  general  conditions  on  the 
ODEs,  the  equilibrium  set  will  be  a  smooth  mathematical  set  called  an  "m-dimensional  mani¬ 
fold."  Although  this  m-dimensional  manifold  is  often  parameterized  by  the  m  control  inputs 
([R],[YSJ],[MKC],(CM],  and  [G])  it  is  sometimes  easier  to  choose  a  different  set  of  coordi¬ 
nates  on  this  manifold.  We  could  choose,  for  example,  some  of  the  state  variables  from  the 
ODE,  or  functions  of  the  state  variables.  Instead  of  first  fixing  the  m  controls  and  calculating 
the  associated  equilibrium  values  of  the  states,  we  can  fix  m  of  the  state  variables  and  solve 
for  the  equilibrium  values  of  the  other  states  and  the  controls. 

For  a  given  set  of  state  variables,  the  dependence  of  the  ODE  on  one  subset  of  the  state 
variables  may  be  more  complicated  than  on  the  others.  For  example,  the  ODE  may  be  tran¬ 
scendental,  discontinuous,  or  even  tabular  in  some  of  the  variables,  while  only  polynomial  or 
even  linear  in  others.  Analytic  computations  will  be  simplest  if  we  choose  as  coordinates  on 
the  equilibrium  manifold  the  m  states  on  which  the  ODE  has  the  most  complicated  depen¬ 
dence.  This  way,  for  each  value  of  these  coordinates,  the  equilibrium  equations  can  be  solved 
by  relatively  simple  means.  If  the  equations  are  algebraic  in  the  remaining  variables,  we  can 
solve  them  using  standard  techniques  [BCLA],[PY],[vdW],[W]  .  In  subsection  4.4,  we  expli¬ 
citly  compute  the  equilibrium  manifold  for  a  six-degree-of-freedom  aircraft  model. 
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4.2  Class  of  Systems 


The  class  of  systems  described  in  this  section  includes  the  aircraft  example  shown  in  Fig¬ 
ure  3 2  which  has  8  states,  4  controls,  and  6  integrator  chains  of  lengths  1,  1,  1,  2,  2,  and  1  . 

Given  a  system  of  ordinary  differential  equations  of  the  form  x  =  F(x,u)  ,  we  can  rewrite 
it  as  F(x^c,u)  =  0  by  setting  F(x,x,u)  =  F(x,u)  -  x  .  We  will  consider  systems  with  n 
states  x,  k  integrator  chains  (degrees  of  freedom),  and  m  controls  u.  The  indexed  set  of  ordi¬ 
nary  differential  equations  can  be  written  as  (see  Figure  4.1) 

^i(Xi^i*X2^j' "  *  ~0  i  —  1,  k  (4.2.1) 

xjj+j  —  Xy  i  —  1, ...,  k  j  —  1, ...,  Cj  —  1 

where  Ci  +  Cj  +  •  •  •  +  ck  =  n 
(  ci  through  ck  are  the  lengths  of  the  k  integrator  chains.) 


x  =  {  xy  }  i  =  1 . k  j  =  l . Cj  for  each 

where  Xy  is  the  j*  integrator  output  (state)  in  the  1th  chain. 


u  =  {  Uj  J  i  =  1,  ...,  m 


We  assume  that  F  has  the  necessary  properties  to  ensure  existence  and  uniqueness  of  solu- 


Typically,  F(x,x,u)  may  be  written  as  F(x,u)  -  x,  and  in  this  case  the  system  may  be 
represented  schematically  as  shown  in  Figure  4. 1 . 


B 


6 
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Integrator  Chains 


Figure  4.1  Schematic  Representation  of  Typical  System 

By  definition,  equilibria  are  solutions  of  equation  (4.2.1)  for  which  x  is  constant.  Conse¬ 
quently  the  equilibrium  set  is  the  set  of  solutions  of  the  following  equations  in  the  (x.u)  space. 

(4.2.2) 

Fi(0,0,  •  •  •  ,0,x,u)  =  0  i  =  I . k 

Xjj*!  =0  i  =  1, ....  k  j  =  1 . Cj  -  1 


This  obviously  reduces  to  k  equations  in  the  m  entries  of  u  and  the  k  nonzero  entries  of  x  , 
(  x, ,  i  =  1,  ....  k  ).  From  the  Implicit  Function  Theorem  [M],  at  the  regular  points,  the 
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equilibrium  space  is  (locally)  an  m-dimcnsional  manifold. 

Sometimes,  by  changing  the  coordinates  of  the  system,  it  is  possible  to  simplify  computa¬ 
tions.  If  F  is  algebraic  in  k  of  the  elements  of  (x,h),  where  (xfx^.u^hfxj^.u))  is  any  change 
of  coordinates  on  the  (x^.u)  space,  then  parameterize  the  equilibrium  set  by  the  other  m 
entries  of  (x(xi>1,u)4i(xi^1,u)).  For  each  value  of  these  m  parameters,  the  equilibrium  computa¬ 
tions  reduce  to  solving  k  algebraic  equations,  F,  in  k  unknowns.  Solutions  can  be  found  by  a 
number  of  methods,  e.g.  repeated  resultants,  Grobner  bases,  etc.  [BCLAJ,[PY],[vdW],  and 
[W].  If  F  is  algebraic  in  more  than  k  of  the  elements  of  (x(x^1,u),h(xi^1,u)),  solve  the  alge¬ 
braic  system  for  the  k  algebraic  elements  of  lowest  degree.  This  reduces  the  number  of  calcu¬ 
lations  involved  in  solving  the  algebraic  system.  If  the  system  is  linear  in  some  of  these  k 
variables,  then  these  can  be  eliminated  by  linear  algebra,  leaving  a  simpler  algebraic  system  to 
solve. 

The  change  of  coordinates  to  (x(xi<1,u),h(xi>1,u))  can  be  essential  in  getting  the  system  into 
algebraic  form,  and  very  helpful  in  reducing  the  order  of  the  system  of  algebraic  equations. 
This  will  be  demonstrated  in  subsection  4.3  . 


4.3  Equations  of  the  Form  F(x,x,u)  =  f(x,x)  +  g(x)  h(x,u) 


Dynamics 

In  the  previous  section  we  were  using  state-space  systems  defined  explicitly  by 

x  =  F(x,u) 


or,  by  letting 


F(x,x,u)  =  F(x,u)  -  x 


the  same  system  was  defined  implicitly  by 

F(x,x,u)  =  0.  (4.3.1) 

This  is  a  system  of  k  equations.  We  would  like  to  analytically  reduce  this  to  a  set  of  k-m 
equations  in  m  fewer  unknowns.  Whether  or  not  we  can  do  this  depends  on  the  algebraic 
structure  of  the  system  of  equations.  In  the  rest  of  this  subsection  we  will  be  considering  a 
case  for  which  this  reduction  is  possible.  Note  that  for  the  aircraft  equations  of  motion  k  =  6 
and  m  =  4,  so  k-m  =  2  . 

Consider  the  case  where  F(x,x,u)  has  the  special  form 

F(x,x,u)  =  f(x,x)  +  g(x)  h(x,u)  (4.3.2) 

Equation  4.3.1  now  takes  the  form 

0  =  f(x,x)  +  g(x)h(x,u)  (4.3.3) 


where 


f  is 

k  x  1 

( column  vector ) 

g 

k  x  m 

(  matrix  with  k  >  m  ) 

h  is 

m  x  1 

(  column  vector  ) 

In  this  case,  the  only  dependence  of  F  on  the  in  controls,  u,  enters  through  the  m  functions 
h(x,x,u).  F  depends  on  h  in  a  linear  way,  so  we  can  eliminate  h  using  linear  algebra.  To 
eliminate  h,  split  the  equations  into  two  parts 


0 

.0. 


Pi  f(x,x)l 

Pi  g0>) 

P2f(x,x)J  + 

P2  SOO 
»  , 

h(x,u) 


(4.3.4) 


where 

Pi  is  any  m  rows  of  a  k  by  k  identity  matrix. 

P2  is  the  remaining  k-m  rows  of  the  k  by  k  identity  matrix. 

If  the  g(x)  matrix  is  full  rank,  it  has  m  independent  rows.  In  this  case  choose  Pt  to  select 
out  these  independent  rows  so  that  Pig(x)  is  invertible,  then 

h(x,x,u)  =  -  (Pig(x))-1  P if(x,x)  (4.3.5) 

Plug  this  expression  for  h  into  the  0  =  P2f(x,x)  +  P2g(x)  h(x,x,u)  equation  to  get 

0  =  P2f(x,x)  -  P2g(x)  (P!g(x))-1  Ptf(x,x)  (4.3.6) 


Separating  out  f  gives 

0  =  (p2  -  P2g(x)  (P ig(x))-1  P,  j  f(x,x) 


(4.3.7) 


or 

0  =  gx(x)  f(x,x)  (4.3.8) 

where  g^(x)  is  the  following  k-m  by  k  matrix  : 

gx(x)  =  (p2  -  P2g(x)  (P^x))-'  P,]  (4.3.9) 

Note  that  (g-L(x)  g(x))  =  0  ,  this  is  the  reason  for  chosing  the  "perpendicular"  superscript  on  g, 
ie  g^.  Equation  4.3.8  can  be  derived  directly  from  equation  4.3.3  by  premultiplying  equation 
4.3.3  with  gJ-(x) . 
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Equilibrium 

At  equilibrium  x  =  0  and  the  only  nonzero  entries  of  x  are  (  Xy  i  =  1,  ...,k  ).  So  at 
equilibrium,  equation  4.3.7  represents  k-m  equations  in  k  unknowns.  Therefore  the  solution 
set  is  m-dimensional. 

Consider  the  case  where  f  is  polynomial  in  k-m  of  the  nonzero  entries  of  x,  and  g  does 
not  depend  on  these  k-m  variables.  Split  the  nonzero  entries  of  x  in  to  two  groups  y  and  z  . 

Let  y  ®  (  *u  }  where  i  takes  on  k-m  of  the  values  that  correspond  to  variables  in  which 
f  is  polynomial. 

Let  z  =  {  Xjj  }  where  i  takes  on  the  remaining  m  values. 

The  function  f  is  polynomial  in  the  entries  of  y.  Let  p  be  the  number  of  monomials  in 
the  entries  of  y  on  which  f  depends,  e.g.  if  f  =  1  +  3yj  +  y2  +  yi7y2  t^icn  P  =  4  ■ 


Let  y  be  a  column  vector  containing  these  p  monomials. 

In  this  case 

f(0,x)  =  f(z)  y  (4.3.10) 

where  f(z)  is  a  k  by  p  matrix  of  coefficients.  See  (4.5.9)  through  (4.5.11)  for  an  aircraft 
example. 

Since  g  only  depends  on  z,  g^  has  the  form 

g^z)  =  [p2  -  P2g(z)  (P1g(z))-'  p,)  (4.3.11) 


and  Equation  4.3.8  becomes 


(4.3.12) 


■  if 


Equation  4.3.12  is  a  set  of  k-m  polynomials  in  the  k-m  entries  of  y,  with  coefficients  that 
depend  on  z  .  For  each  value  of  z,  this  set  of  polynomials  can  be  solved  using  repeated  resul¬ 
tants  or  various  other  methods  [BCLA],  [PY],  fvdW),  and  (W).  These  techniques  reduce  the 
system  of  k-m  polynomials  in  k-m  variables  into  a  new  system  of  k-m  polynomials  of  the  fol¬ 
lowing  form. 


i 

i 

sa 
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<ik-m(yi.y*  •  *  •  .yk-m)  =  0 


qk-m-i(yi.y2.  *  *  ov-m-i)  —  ® 


<h(y\*yi)  =  o 

qi<yi)  =  o 

The  qi(yi)  polynomial  in  one  variable  can  be  solved  numerically  by  placing  its  coefficients 
into  a  companion  matrix,  then  finding  the  eigenvalues  of  this  companion  matrix  using 
EISPACK  [SBDGIKM],  If  the  qjCyj)  polynomial  is  represented  by 


qi(yi)  -  *0  +  »iyi  +  •  •  •  +  +  w* 

then  the  companion  matrix  is  of  the  form 


(4.3.13) 


0 

1 

0 

...  0 

0 

0 

1 

...  0 

0 

0 

0 

... 

• 

• 

• 

...  0 

0 

0 

0 

...  1 

_  i2. 

_  ai_ 

_  ^2 

_  a^~l 

ad 

«d 

ad 

(4.3.14) 


The  eigenvalues  of  this  matrix  will  be  the  roots  of  the  polynomial  q{  .  Each  root,  y1  ,  of  this 
polynomial  is  then  plugged  back  into  the  q2(yi»y2)  polynomial.  Since  y,  is  now  known, 
q2(yi.y2)  *s  now  a  polynomial  in  only  one  unknown,  y2  .  Put  the  coefficients  of  the  q2  poly¬ 
nomial  into  a  companion  matrix  and  solve  for  y2  .  Continue  this  process  of  back  substitution 
until  all  the  values  of  (yj,y2,  •  ■  •  ,y k_m)  are  solved  for. 

At  this  point,  y  has  been  determined  for  a  given  set  of  values  of  z.  Since  the  nonzero 
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components  of  x  consist  of  the  entries  of  y  and  z,  we  now  have  all  the  entries  in  x,  so  we  can 
solve  for  the  inputs,  u,  by  solving  equation  4.3.5 

h(x,u)  =  -  (P1g(x))-1  P1f(0,x) 

for  u  . 

All  the  points  in  the  equilibrium  space  can  be  calculated  by  dividing  the  m-dimensional 
set  of  coordinates,  z,  into  a  grid.  For  each  point  in  this  grid,  calculate  the  corresponding  y 
and  u.  An  (m- 1  )-diinensional  family  of  one-dimensional  slices  of  this  equilibrium  space  can 
be  plotted  to  visualize  the  space. 

The  methods  for  solving  systems  of  polynomial  equations  work  in  principle  for  any  number  of 
polynomials,  they  are  much  simpler,  however,  when  the  number  of  polynomials  is  small.  For 
example,  a  system  of  two  or  three  polynomials  each  of  degree  two  or  three  can  be  solved 
quickly,  while  a  system  of  four  or  five  polynomials  each  of  degree  two  or  three  can  take 
several  minutes  of  computer  time  to  solve. 

So  far  the  methods  described  are  quite  general.  In  the  next  section,  we  will  be  consider¬ 
ing  a  concrete  example.  For  the  aircraft  example,  k  =  6  and  m  *  4,  so  we  are  left  with  2 
equations  in  6  unknowns.  The  equations  are  polynomial  in  several  of  the  unknowns.  Choose 
the  two  (k-m)  entries  of  y  from  these  polynomial  type  unknowns.  Use  the  remaining  four 
unknowns  (  z  )  to  parameterize  the  equilibrium  manifold.  For  each  set  of  values  of  these  four 
parameters,  solve  the  two  polynomial  equations  for  the  remaining  two  unknowns  (  y  )  . 


4.4  Introduction  to  Aircraft  Equilibrium 


The  kind  of  equilibrium  manifold  that  we  get  will  depend  on  what  we  choose  for  the  states 
(outputs  of  integrators).  We  will  not  include  ?  as  a  state  since  none  of  the  forces  or 
moments  depend  on  4*  .  At  equilibrium,  4*  can  take  on  any  constant  value  (including  zero), 
so  the  equilibrium  trajectories  will  be  vertical  helices  (see  Appendix  A).  Note  that  straight 
line  trajectories  are  just  infinitely  fat  vertical  helices  with  4*  =  0  .  "Steady"  maneuvers  such 
as  a  3  g  pullup  arc  precluded  since  the  component  of  gravity  is  varying  greatly  during  such  a 
maneuver.  It  would  still  be  possible  to  trim  about  a  3  g  pull  up  if  we  considered  time  varying 
ODE  and  trimmed  about  some  nominal  trajectory.  When  we  refer  to  equilibrium  in  this 
report,  we  will  be  referring  to  the  vertical  helix  type  of  equilibrium. 

Equilibrium  is  attained  when  the  six  outputs  (V,(i,a,<D,0,4r)  all  remain  constant  .  If  we 
specify  the  constant  values  for  four  of  the  outputs  then  the  equilibrium  equations  will  deter¬ 
mine  the  corresponding  values  for  the  four  inputs  and  the  remaining  two  outputs. 

For  several  choices  of  sets  of  four  constant  outputs,  the  equations  are  algebraic  in  the 
remaining  variables.  This  allows  us  to  solve  the  system  of  nonlinear  equations  for  the 
corresponding  inputs  and  the  remaining  two  outputs.  Two  cases  will  be  considered. 


Case  1:  The  simplest  case  occurs  for  low  speed  flight  (say  Mach  <  .6)  where  the  aero¬ 
dynamic  functions  are  independent  of  Mach.  For  instance  Cx(M,a,p)  becomes  just  Cx(a,(3)  . 
In  this  case  the  equations  are  second-order  polynomial  in  V  and  'F  ,  so  it  is  easiest  to  specify 
constant  values  for  the  four  outputs  (|3,a,<X>,9)  then  solve  the  system  of  two  second-order 


polynomials  for  the  remaining  two  outputs  (V,40  and  then  solve  for  the  inputs 


T 

5. 

5. 


Case  2:  For  high-speed  flight,  when  the  aerodynamic  functions  depend  on  Mach  number 
(Mach  is  a  function  of  V  and  air  temperature)  the  equations  are  no  longer  polynomial  in  V  . 
However,  the  equations  are  still  second-order  polynomial  in  4*  and  second-order  polynomials 
in  trigonometric  functions  of  d>  and  0.  Since  trig  functions  are  rational  functions  of  complex 
exponentials,  the  equations  are  rational  in  e1<1>  and  e1®.  By  multiplying  the  equations  through 
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by  the  common  denominator  of  these  rational  functions,  the  equations  become  fourth-order 
polynomial  in  e*  and  e*®.  In  this  case  it  is  easiest  to  specify  constant  values  for  the  four  out¬ 
puts  (V,p,a,<I>)  then  solve  a  system  of  two  polynomials  for  the  remaining  two  outputs 

Tt" 


(e*®,'F)  and  then  solve  for  the  inputs 


8. 
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4.5  Calculation  of  Aircraft  Equilibrium 


At  equilibrium,  x  =  constant;  Le.  (V,p,a,<j>,0,H/,<I>,0)  =  constant  Since  0,0  are  constant, 
0  =  0  and  <i>  =  0  ,  so  (3.3.2)  becomes 


pi  1  0  -sin(0)  <j> 

Q  =  0  cos(<D)  cos(0)sin(<b)  0  = 

0  -sin(<D)  cos(0)cos(<l>)  'P 


(4.5.1) 


1  0  -sin(0) 

0  cos(G>)  cos(0)sin(<t>) 
0  -sin(d>)  cos(0)cos(<I>) 


o]  .  r 

0  =VF  c 
V  IP 


-sin(0) 

cos(0)sin(d>) 

pos(0)cos(<X>) 


This  says  that  JQ  J  is  constant  so  the  —  operator  in  the  covariant  derivative  of  the  angular 

[VI  fu” 

momentum  can  be  eliminated.  Similarly,  p  is  constant  which  implies  V  is  constant  so 

laj  LW. 

the  operator  in  the  covariant  derivative  of  the  linear  momentum  can  be  eliminated.  Also 

U  cos(P)cos(a) 

note  that  in  general  (3.3.3)  gives  us  V  =  V  sin(p) 

[w  lcos(p)sin(a) . 

Using  the  expressions  for  P,  Q,  and  R  from  equation  4.5.1,  the  skew  symmetric  matrix  in 
the  covariant  derivative  expression  becomes: 


0  R 

-Q 

[  0 

cos(0)cos(d>) 

-cos(0)sin(d>) 

-R  0 

p 

=  *P  -cos(0)cos(<X>) 

0 

-sin(0) 

.Q  -P 

0  . 

cos(0)sin(<I>) 

sin(0) 

0  J 

We  will  use  the  following  names  to  keep  the  notation  more  compact: 


1 

$ 

jC 

[rt 


B 


I 


cos(0)cos(<I>)  -cos(6)sin(<I>) 


E*e  =  -cos(0)cos(<I>) 

0 

-sin(0) 

[  cos(0)sin(O) 

sin(0) 

0 

fcos(P)cos(a) 

[  -sin(0) 

lpa  =  sin(p) 

^4>e  = 

cos(0)sin(<I>) 

Lcos(p)sin(a) . 

lpos(0)cos(<t>) 

Qyi  “ 


On.  Q*<}  Qi, 

=  Cy<3  ^y* 

pT  OtQ 

I„  -I,y  -IU 

^mo  =  “ ^*y  Vy  — 

~^X2  — ^yz  ^Z2 


Qmn  " 


Cl»  C1q  Cl* 


CinrnfQ*  —  On*  ^m<3 

Q|>  ^"q 


b  0  0 
BCB  =  0  c  0 
0  0b 


V2  _  mg 

nom"  VipS 


Plugging  all  this  into  equations  (3.3.1),  (3.3.5),  (3.3.7),  (3.3.10),  and  (3.3.4)  gives  the 
equilibrium  expressions  for  f,  g,  and  h  in  terms  of  (V,p,cx,'I/,<J>,0)  and  the  four  controls. 

The  top  three  rows  of  f(0,x)  become 

W 

-j-E^lpa  +  W  + 


y  ^*yz  +  2V  ^y***®^  1<1>e 

v  nom  ^  T 


(4.5.2) 


while  the  bottom  three  rows  of  f(0,x)  become 


g 

i 


(4.5.3) 


~rBCB'lEWmo1«e  + 


¥ 


^lran  "*■  2V  Qntn»q*B^B  ^♦O 


g(x)  is  left  unchanged  since  each  entry  only  depends  on  (M(V),p,a) 

^*t  ^  C**'  C*t, 


g(x)  = 


S  Cy^  Cyt 


<*  cu 


ck  <v 


Qi»t  Qn«t 

Qvr  ^»U  ^n«  ^n* 


and  the  h(0,x,u)  vector  becomes 


h(0,x,u)  = 


,  m8  . 

V2  .  Vsm(e)b. 
_  sin(5a+ - — — ) 

*  nnm  ** 


V2 


V2„om 

V2 


V2, 


sin(8e-a) 

sin(5r-0) 


nom 


(4.5.4) 


The  g^(x)  matrix  described  in  equation  (4.3.9)  requires  the  selection  of  Pj  and  P2  .  The 
choice  given  below  is  motivated  by  the  desire  to  make  [Pj  g(x)]  invertible  for  typical  aircraft. 
Let 


P.  = 


1  0  0  0  0  0 
0  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 
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Let 


0  1  0  0  0  0 
0  0  1  0  0  0 


Using  Plt  P2,  and  g(x)  we  can  form 

gHx)  =  P2  -  PigCxKP^Cx))"1?, 


(4.5.5) 


Note:  For  aircraft  with  left/right  symmetry,  the  elevator  and  thrust  only  affect  the  first,  third, 
and  fifth  rows  of  f(x,x)  +  g(x)  h(x,x,u)  while  the  rudder  and  aileron  only  affect  the  second, 
forth,  and  sixth  rows  of  f(x,x)  +  g(x)  h(x,x,u)  .  The  result  is  that  the  g(x)  matrix  then  has  the 
following  simplified  form: 


g(x)  = 


CxT 

0 

0 

0 

cy*. 

0 

°y‘, 

CIx 

0 

c*. 

0 

0 

0 

cu. 

CmT 

0 

Cn* 

0 

0 

Cn*. 

0 

s 

(4.5.6) 


In  this  simpler  situation,  P^x)  can  be  inverted  explicitly  so  equation  (4.5.5)  gives  gJ*(x)  as  : 

Cy,  Cn,  -  Cy,  Cy,  Ci.  “  Cy,  Q. 

0  1  0  - - — - : — -  0  - - — - - -1 

C^Cn.-  Cl,  ^ 

=  Cz.  C^-  Cz.C*.  Cx,  Cz,-  CXtCz, 


CXtQii*  “  CX«  Qnr 


Cx.  Cz,-  CXtCz^ 
CxTCjn.  -  CXj  CnT 


In  the  further  simplifying  situation  where  the  thrust  is  aligned  with  the  x  axis  (C7t  =  0  =  C^) 
and  the  ailerons  only  produce  a  rolling  moment  (Cy6  =  0  =  ),  we  get: 


0  10  0  0 


gx(x)  = 


0  0  10- 


Cy«. 


(4.5.8) 


If  Cy4  and  Cu  arc  also  neglected,  then  =  P2  .  In  general  the  g(x)  matrix  can  be  a  full 

matrix,  and  g^(x)  can  be  computed  numerically  for  each  value  of  x  used.  In  the  case  where 
g-*-(x)  is  computed  numerically,  it  is  better  to  use  the  singular  value  decomposition  algorithm 
(e.g.  in  UNPACK  [DMBS]  )  in  place  of  equation  (4.5.5)  because  the  SVD  algorithm  is 
numerically  more  stable  and  works  even  when  Ptg(x)  is  not  invertible. 


We  now  have  two  equations, 

-  g1(x)f(0,x)  (4.5.9) 


in  the  six  unknowns,  x  =  (V,  p,  a,  41,  9,  <X>) . 


Case  1  :  (same  as  in  subsection  4.4)  The  case  where  the  aero  functions  do  not  depend  on 
Mach. 


Let  y  =  (V.'P),  y  = 


1 

VV 

V'P 
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,  and  z  =  (P,a,<D,9)  . 


In  this  case  f(z)  is  the  following  6  by  4  matrix,  where  O3  is  a  column  of  three  zeros. 

C 


f(z)  = 


<te 


-xyz 


c*yzfQ«  BCB  .  E<t>e  . 

+ - 1 


V* 


0, 


'Imn 


V2 


^Ijnnpq^  BCB 
V2  2 

nom 


g 


l<t>e 


BCB~ 

mg 


O3 


EWmol<&e 


so  the  equilibrium  equations  become 


33 


=  g^p.a.G.ewp.a,*,©)  ^ 


Choose  values  for  P,a,<X>,6  and  solve  the  above  set  of  two  equations  for  V.Y.  These  equa¬ 
tions  are  second -order  polynomials  in  V  and  4*  ,  of  the  following  form: 


0  =  pOVb  =  Poo  +  Pjo  Vj  +  p„  V  ¥  +  Poj  <i-2 


0  -  q(V,<h  =  <ta)  +  q»  V2  +  q„  V  >}'  +  qoj  >P3 


(4.5.10) 


(4.5.11) 


where 


Poo  P20  P11  P02 
Q20  Qll  ^02 


=  ^(p.a,*.©)  f(p,a,d>,e) 


The  two  polynomial  equations,  (4.5.10)  and  (4.5.11),  can  be  solved  using  direct  elimination, 
or  by  using  resultants.  Two  polynomial  equations  (whose  leading  coefficients  are  nonzero) 
have  a  solution  if  and  only  if  their  resultant  is  zero  [vdW]  .  Given  two  polynomials  polyl 
and  poly2  in  several  variables,  of  degree  nl  and  n2  respectively  in  some  particular  variable 
(say  x),  their  resultant  with  respect  to  x  is  the  determinant  of  the  (nl+n2)  by  (nl+n2)  matrix 
formed  as  follows.  Take  the  coefficients  on  powers  of  x  from  polyl  and  put  them  in  the  top 
row  of  the  matrix  starting  on  the  left  side,  repeat  this  row  n2  times,  shifting  to  the  right  by 
one  column  for  each  row.  Fill  in  row  n2  +  1  with  the  coefficient  on  the  powers  of  x  from 
poly 2.  Repeat  this  row  nl  times,  shifting  to  the  right  by  one  column  for  each  row. 

The  resultant  of  polynomials  (4.5.10)  and  (4.5.11)  (with  respect  to  4*  with  Po2qo2  *  0  )  is 
the  determinant  of  the  following  matrix  whose  entries  are  shifted  rows  of  the  coefficients  on 
the  4*  terms.  (  nl  +  n2  =  4  unless  p02  or  q02  is  zero). 


i 

fcimiranrr  rrr  -rnmn-Tr 


ViV.V.V.V.%  .V 


g 


(4.5.15) 


0  =  dct 


Poo+P^V2 

0 

qOO-K^oV2 

0 


PllV  P02  0 

POO+PttV2  Pliv  P02 

qnV  qo2  o 

qoo+<hoV2  qyV 


k4V4  +  k2V2  +  ko 


k4  =  (P2oqo2"q2oPo2)2  ~  (Pnqo2-qnPo2)(P2oqn^q2oPti)  (4.5.1 6) 

^2  =  2  (Pooq02^00P02XP20q02"^20P02)  “  (Pliq02~qtlP02)(Pooqil~qooPll)  (4.5.17) 

ko  =  (Pooqo2~qooPo2)2  (4.5.18) 

The  k’s  depend  on  the  p^'s  and  q^'s  so  they  depend  on  (5,a,<I>,0  .  Since  no  odd  order 
terms  in  V  appear  in  the  polynomial  on  the  right  hand  side  of  equation  (4.5.15)  ,  it  is  of 
second-order  in  V2  so  it  can  be  solved  using  the  quadratic  formula: 


2  -k2  ±  V(k2k2  -  4k4k0) 

2k4 


(4.5.19) 


From  (4.5.19)  we  see  that  there  are  at  most  two  positive  real  solutions  for  V.  From 
numerical  evaluation  of  several  aircraft  models,  we  have  seen  that  typically  one  solution  of 
(4.5.19)  is  large  and  positive  while  the  other  is  so  small  that  the  aerodynamic  control  surfaces 
would  not  have  enough  authority  to  hold  the  aircraft  at  that  equilibrium.  The  equilibrium 
associated  with  the  large  speed  is  typically  associated  with  low  'V  values  and  corresponds  to 
reasonable  flight  regimes,  while  the  low-speed  equilibrium  is  often  associated  with  large-4/, 
spin-type  flight 

Use  the  real  positive  values  of  V  to  find  ¥  as  follows.  First,  eliminate  4/2  from  the  two 
polynomials  by  forming  the  expression 

q02  P^)  -  P02  q(V,40  • 

The  S'2  terms  cancel  in  this  expression,  leaving 


0  =  qo2(Poo+P20v2+Pi  i ~  P(tt(<toorKl20^2+<lii'r'*/) 


(4.5.20) 


,j,  _  (q02POQ~P02qoo)  +  (q02P2Q~P02<ho)v2 
(P02<lll“<l02Pll)V 


(Note:  If  qo2  and  are  both  zero,  then  there  is  no  'F2  term  to  stan  with.) 
Use  these  values  of  V  and  *F  in  f  and  g  to  solve  for  h  : 


(4.5.21) 


h(0,x,u)  =  4P1g(x)]~1  [P,f(0,x)J 


g 

Finally,  use  this  value  of  h  to  solve  the  following  equations  for  g 


v2  .  /C .  ^si^ejb 

__  s,n<8.+  ' 

^  nom 

V2 

— —  sin(8e-a) 

'  nom 

V2 

~2 —  sin(8r-P) 


=  h(0,x,u) 


Case  2  :  (See  subsection  4.4)  When  the  aero  functions  depend  on  Mach  number,  let 
z  =  (V,  P,  a,  <I>) ,  y  =  (el0,'F)  ,  and  $  is  the  vector  of  7  monomials  in  e,e  and  'F  on  which  f 
depends: 

9  =  transpose(  (e_10,l,e,0),H/(e_ie,eie),vi/2(e_i2e,l,ei2e)  )  . 


The  two  equations  (4.5.9)  are  second-order  polynomials  in  4*  ,  sin(@)  ,  and  cos(©)  .  The 


v; 

Sf 

5*! 

S 

ft 

ft 


*•.< 

s< 


y. 

I**' 

wl 

I 

& 

1 

•M 

$ 

K 


SI 

a 

SI 

$ 

il 

* 


V*.m*G*  W  JV \TV 


only  0  dependence  comes  from  the  and  in  f(0,x)  (see  equations  4.5.2  and  4.5.3) 
Substitute  for  cos(6)  and  sin(0)  using 


sin(0)  = 


eie  -  e~*9 


cos(0)  = 


eie  +  e“ie 


where  i  =  V(-l)  . 


i®  ../i  --i® 


This  lets  us  write  and  E$e  in  terms  of  e1**  and  e 


=  E^  ei0  +  E^  e 


1<^0  —  1®  eie  +  1<£  e**0 


where 


n  cos(<D)  -sin(<I>) 
0  2  2 

,-cos.^l  0  zL 

^2  2i 

sin(<D)  J_  Q 

2  2i 

and  the  over-bar  signifies  complex  conjugation. 


♦  1  $ 


-1 

2i 

_  sin(<I>) 
2 

cos(<I>) 

2 


Plug  these  expressions  for  1M  and  into  (4.5.2)  and  (4.5.3)  .  The  result  is 


where 


f(0,x)  =  f(z)  y  =  ff123(z)  .  f45(z)  .  f678(z>]  ^*45 

***678 


9m  =  I 


Y45  -  ei0 


*678  =  1 


v*»*. .'I  -•»**•* Vi  A*l*k*i *.***.*«* 


^123(z)  = 


■CTV,  1 


tyz 


°3  ^2 - Clmn  °3 


•~E«1Pa  + 


V  r  BCB  t 
v2  *y>o«  2  * 


^4S(Z)  = 


V  V 

+  "y2 


V2  2 
T  nom 


V  r  BCB  - 

V2  ^hwVQ,  j 


r  BCB, 

Minn**  j  ** 


^67«(Z)  = 


mg 


^3 

mg 


BCB-1 

mg 


Multiply  both  sides  of  the  equations,  O2  *  g^zjftz^  ,  by  ea®  to  make  y  ei20  polynomial  in 
e*®  .  The  resulting  two  equations  are  second-order  polynomial  in  and  fourth-order  polyno¬ 
mial  in  e1®  of  the  form  : 


0  =  pote*®)  +  p1(eie)'F  +  p2V9)'V2 


(4.5.22) 


0  =  q0(O  +  q1(e,w)^/  +  q2(elW)'P 


0\vi/2 


(4.5.23) 


where 


2^  *  ^(Z)  f|23(z)  [yi23  Ci20l 


q,1^®)  =  gl(z)  "f4s(z)  ^4S  ^ 


[• 

& 

ijj'i 

ft* 

1 


I 
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J 

5 

9 

«» 

6 

!$ 

$ 

J 


ft 


1 

I 


m 


II 


Wit 


H 


2S 


Form  the  resultant  (with  respect  to  4*)  again 

Po  Pi  Pi  0 

0  Po  Pi  Pi  , 

0  =  dCt  qo  qi  qi  0  =  (Po^-P^o)  -  (Piq2-P2qi>(Pofll-Piqo)  • 

0  <k>  4i  42 

The  right  hand  side  will  be  a  14th-order  polynomial  in  e,e  of  the  form  : 

e“®  £  V**  (4.5.24) 

fc— 6 

where  c_k  =  ck  (so  wc  only  need  to  calculate  Cq  through  c6  ).  The  polynomial  in  (4.5.24)  is  of 
the  form  (e*2®)  (12th-onder  polynomial  in  e*®  )  .  Put  the  coefficients  of  this  12th -order  poly¬ 
nomial  into  a  companion  matrix  and  find  the  roots  using  EISPACK  [SBDGIKM].  Keep  the 
roots  that  have  magnitude  equal  to  1  (6  is  real  in  these  cases).  Use  these  roots  to  solve  for 
4*. 

po(ei0)  ~  P2<eie)  qo(eie) 

P2(cie)  qi(eie)  -  q2(e*®)  Pi(eie) 

4*  will  be  real  since  the  extra  factors  of  e,e  cancel  in  numerator  and  denominator. 


Finally,  solve  for  the  inputs  as  in  case  1  . 


'.wvwvunwwjw 


SECTION  5  :  ANALYSIS  OF  AIRCRAFT  DYNAMICS 

In  this  section  we  analyze  the  nonlinear  aircraft  equations  using  dynamic  inversion.  Our 
goal  is  to  show  how  dynamic  inversion  can  be  used  to  develop  flying  quality  parameters. 
Some  candidate  parameters  arising  from  the  analysis  in  this  section  are  presented  in  section  6: 
Nonlinear  Flying  Qualities. 


a 
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5.1  Partial  Dynamic  Inversion 


The  6  DOF  nonlinear  aircraft  equations  of  motion  are  of  the  form: 


x  =  f(x)  +  g(x)  h(x,u)  (5.1.1) 

with 

x  in  an  n  dimensional  state  space  (n»S  for  our  aircraft  model) 
u  in  an  m  dimensional  control  space  (m=4  for  our  aircraft  model) 
f  :  an  n  by  1  vector  depending  on  x 
g  :  an  n  by  m  matrix  depending  on  x 
h  :  an  m  by  1  vector  depending  on  x  and  u 


Equation  5.1.1 


splits  into 


the  following  form 


•  ■ 

*1 

fj(x) 

gl(x) 

*2 

s 

f2(x) 

•f 

g2(x) 

*3 

•  * 

f3(x) 

0 

h(x,u) 


(5.1.2) 
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P 


•  .*  «l)  < .«  t.*  nt  v.a  t.|  a.a  >  *  j  , 


I  itl'M  t 


where  the  first  k  equations  represent  the  dynamics  (k=6  for  a  6 DOF  aircraft): 

*1  fi(x)  gi(x) 

-  =  -  +  -  h(x,u) 

x2  fjOO  g2(x) 


(5.1.3) 


and  the  last  n-k  equations  represent  the  kinematics  (n-k=2  for  our  6DOF  aircraft  model): 


x3  =  f3(x) 


(5.1.4) 


We  split  the  dynamic  state  in  equation  5.1.3  into  the  xx  and  x2  parts  to  distinguish  the 
states  xt  we  wish  to  control  directly,  from  the  states  x2  we  choose  not  to  control  directly. 
When  m  =  k  we  can  control  all  k  of  the  degrees  of  freedom.  When  m  <  k  ,  we  can  only 
control  pan  of  the  dynamics,  because  there  are  fewer  independent  degrees  of  control  authority 
than  degrees  of  freedom  .  In  equation  5.1.2  we  have  : 

Xj  in  an  m-dimensional  space 
x2  in  a  (k-m)-dimensional  space  (  k  £  m  ) 
x3  in  an  (n-k)-dimensional  space  (  n  £  k  ) 
ft  :  an  m  by  1  vector  depending  on  x 
f2  :  a  k-m  by  1  vector  depending  on  x 
f3  :  an  n-k  by  1  vector  depending  on  x 
g!  :  an  m  by  m  matrix  depending  on  x 
g2  :  a  k-m  by  m  matrix  depending  on  x 


In  cases  where  we  cannot  invert  all  of  the  dynamics,  we  can  still  do  a  partial  inverse. 
The  xt  dynamics  can  be  inverted  as  follows.  Let  v  be  some  function  of  the  pilot’s  com¬ 
mands,  then  given  a  set  of  desired  dynamics  xt  =  Ft(x,v)  ,  put 

h(x,u)  =  (gjtoJ-'t  xt  -  f,(x)  1  =  (g^x)]"1!  Fj(x,v)  -  ft(x)  ]  (5.1.5) 

and  then  solve  h(x,u)  for  u,  the  signal  to  the  actuators,  as  a  function  of  x  (which  we  assume 
can  be  measured)  and  the  external  input  v  .  Equation  5.1.5  involves  the  inverse  of  g[  (x) 
which  might  only  exist  in  some  subset  of  the  state  space.  In  this  sense,  the  analysis  is  not 
global. 


41 


To  determine  the  uncontrolled  x2  dynamics  and  the  kinematics,  substitute  the  expression 
for  h(x,u)  from  5.1.5  into  5.1.2  to  obtain 


x,  =  Fjfov) 

x2  =  f2(x)  +  g2(x)[gi(x)]~1[  F^x.v)  -  f^x)  ] 
x3  =  f3(x) 


This  partial  inversion  process  involves  several  choices.  The  first  choice  made  (choice  of 
xt)  determines  which  dynamic  states  are  controlled  directly.  The  second  choice  made  ( 
choice  of  F((x,v)  )  determines  the  dynamics  for  those  states.  Some  examples  we  explored  on 
this  program  are  shown  below. 


Example  1  :  The  U,P,Q,R  Inverter 


[U 


Fi(x.v)  = 


*2=  [w] 


I*] 

x*=  LeJ 


u  -  uOTd‘ 

P-P«nd 
Q  ~  Qand 
„R  ~  Rcmd . 


IRJ 

Xu  0  0  0 

0  Xp  0  0 
0  0  Xq  0 
0  0  0  XR 


where 


Example  2  :  The  U,P,<I>,0  Inverter 


where 


F,(x.v)  = 


where 


_  a  _  fip] 
x2  =  \j/  x3  ~  [q  J 


Fi(x,v)  = 


Xjj  0  0  0  0  0  U  “  ^emd 

0  \*  0  0  0  0  ^  ~^cmd 

<l> 

0  0  —£4(0$  0  -<d^  0  0 

0  0  0  “Ce<Be  0  -g>4  0  _  q** 


Example  3  :  The  P,a,C>,0  Inverter 


y  _  [ol 

X2  x3  —  [iQj 


Xp  0  0  0  0  o]P-PcnKl 

o  x*  o  000  a"^tcmd 

0  0  0  -gj^  0  0 

000  -<90)9  0  -<4  Zzt'md 


i 
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Example  4  :  An  alternative  P,a,<J>,0  controller 


Given  an  external  command 

Pen*' 

“and 
V"  *cmd 
®cmd. 


find  the  associated  equilibrium  values  of  . 

(W  • 

where  a,  b,  and  c  are  functions  of  pcmd  ,  acmd  ,  ,  9emd  (see  equation  4.5.19  in  the  equili¬ 

brium  section)  . 

If  there  are  two  positive  real  equilibrium  speeds  to  choose  from,  the  larger  one  is  chosen. 
Next  is  computed  (  cosfa^^osCp,.^)  )  .  Finally,  this  Ucmd  ,  along  with 

Pcrnd-  ^cmd-  40(1  •  ecmd  sent  to  the  U,p,d>,0  inverter. 

Example  5  :  A  P, 0,741  controller 

Given  an  external  command 

Pcmd 

_  ®cmd 
V=  Ycmd 
M’cmd . 

where  7  is  the  velocity  pitch  angle  (flight  path  angle)  and  p  is  the  velocity  roll  angle  (rol¬ 
ling  around  the  velocity  vector).  The  P, 0,741  controller  first  performs  a  change  of  coordi¬ 
nates  from  P,^  ,  ocmd  ,  7cmd  ,  p^  to  P^  ,  ^  ,  0^  ,  using  the  following 


transformations: 


I  |.t  I,*  JJJ.*  M  fe* :  L*  .1 


b_sin(0)  cos(a)cos(P)  -cos(a)sin(p)  -sin(a)  -sin(y) 

;(©)sin(<I>)  =  sin(P)  cos(P)  0  cos(Y)sin(ji) 

(0) 


n(a)  -sin(Y) 

0  cos(Y)sin(ji) 

e/n\  bos(Y)COSOl). 


lpos(0)cos(<D) J  |sin(a)cos(P)  -sin(a)sin(P)  cos(a)  J  ^os^^cos^ ■* 

The  top  row  of  this  equation  gives  us 

0  =  sin-1  |cos(a)cos(P)sin(Y)  +  cos(a)sin(P)cos(Y)sinOi)  +  sin(a)cos(Y)cos(ji)  | 
while  the  ratio  of  the  second  and  third  tows  gives  us 


d>  =  tan-1 


)sin(Y)  +  cos(p)cos(Y)sin(u)  + 


-sin(a)cos(P)sin(Y)  +  sin(a)sin(P)cos(Y)sin(ji.)  +  cos(a)cos(Y)cos(p) 


These  values  of  Pcmd  ,  ,  <Dcmd  ,  0^  are  then  used  to  find  the  associated  equili¬ 

brium  values  of  Vend  . 

rv  \2  _  -b±Vb^-4ac 
(Vcmd  _  2a 


where  a,  b,  and  c  are  functions  of  P^  ,  <t>, 

brium  section)  . 


‘cmd  •  ©and  (see  equation  4.5.19  in  the  equili- 


If  there  are  two  positive  real  equilibrium  speeds  to  choose  from,  the  larger  one  is  chosen. 
Next  Ucmd  is  computed  (  cos(Pand)cos(acm<j)  )  .  Finally,  this  Ucmd  ,  along  with 

Pcmd*  ^cmd-  31,(1  •  ecmd  are  scnt  t0  thc  U,p,O,0  inverter  . 


For  an  introduction  to  the  partial  inversion  method,  and  its  application  to  the  3DOF  long¬ 
itudinal  axis  of  an  aircraft,  see  (El]  .  A  related  method  of  controlling  nonlinear  systems  can 
be  found  in  [HSM]  . 

When  using  dynamic  inversion,  the  Xj  vector  contains  m  states  which  can  be  controlled 
exactly  by  the  m  control  inputs.  We  have  no  direct  control  over  x2  and  x3  ,  so  we  have  to 
check  that  they  remain  stable.  This  analysis  is  done  in  subsections  5.2  through  5.5. 


m 
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5.2  Theory  of  Complementary  Dynamics 


We  will  be  considering  systems  with  m  controls  and  k  chains  of  integrators  (with  k  >  m). 
Note  that  for  conventional  6DOF  aircraft  models,  k=6  and  m  =  4  .  When  the  m  controls  are 
used  to  control  m  of  the  integrator  chains,  the  dynamics  of  the  remaining  k-m  chains  will  be 
completely  determined.  We  will  derive  the  expression  for  the  dynamics  of  these  remaining 
integrator  chains. 

We  are  working  with  state-space  systems  defined  implicitly  by  equation  2.1  which  con¬ 
tained  the  relation 

F(x,x,u)  =  0.  (5.2.1) 

This  is  a  system  of  k  equations.  We  would  like  to  analytically  reduce  this  to  a  set  of  k-m 
equations  in  m  fewer  unknowns.  Whether  or  not  we  can  do  this  depends  on  the  algebraic 
structure  of  the  system  of  equations.  In  the  rest  of  this  subsection  we  will  be  considering  a 
case  for  which  this  reduction  is  possible. 

Consider  the  case  where  F(x,x,u)  has  the  special  form 

F(x,x,u)  =  -  e(x)x  +  f(x)  +  g(x)  h(x,u)  .  (5.2.2) 

Equation  5.2.1  now  takes  the  form 

c(x)x  =  f(x)  +  g(x)h(x,u)  (5.2.3) 


where 


e  is 

k  by  k 

f  is 

k  by  1 

g  is 

k  by  m 

h  is 

m  by  1 

(  square  matrix  ) 

(  column  vector ) 

(  matrix  with  k  >  m  ) 
( column  vector ) 


In  this  case,  the  only  dependence  of  F  on  the  m  controls,  u,  enters  through  the  m  functions 
h(x,u).  F  depends  on  h  in  a  linear  way,  so  we  can  eliminate  h  using  linear  algebra.  To  elim¬ 
inate  h,  split  the  equations  into  two  parts 


’j*~r  V-Vr*nww,TH''^,"^;vu#WTrJir^»v."^ 


Pj  C(x)x  P,  f(x)  P,g(x) 

P2  c(x)x  =  P2f(x)  +  P2g(x)  h(x>u) 


(5.2.4) 


where 


Pt  is  any  m  rows  of  a  k  by  k  identity  matrix. 

P2  is  the  remaining  k-m  rows  of  the  k  by  k  identity  matrix. 

If  the  g(x)  matrix  is  full  rank,  it  has  m  independent  rows.  In  this  case  choose  P,  so  that 
Pjg(x)  is  invertible,  then 

h(x,u)  =  (Pig(x))-1  Pt(e(x)x  -  f(x))  .  (5.2.5) 

Plug  this  expression  for  h  into  the  P2e(x)x  =  P2f(x)  +  P2 g(x)  h(x,u)  equation  to  get 

P2e(x)x  =  P2f(x)  +  P 2g(x)  (Pig(x))-1  P,(f(x)  -  e(x)x)  .  (5.2.6) 


Separating  out  e  and  f  gives 


[p2  -  P2g(x)  (Pjg(x))-1  Pt]  e(x)x  =  [p2  -  P2g(x)  (Plg(x))-1  Pt]  f(x) 


.  (5.2.7) 


(5.2.9) 


g^x)  e(x)x  =  g’^'(x)  f(x)  (5.2.8) 

where  is  the  following  k-m  by  k  matrix  : 

g1(x)  =  [p2  -  P2g(x)  (Plg(x)rl  P,  ]  (5.2.9) 

Note  that  [g^MgCx)!  =  0  ,  so  equation  5.2.8  can  be  derived  directly  from  equation  5.2.3  by 
premultiplying  equation  5.2.3  with  g-L(x)  . 

x  is  the  column  vector  with  elements  (Xjci)  i  =  l,  ....  k  (see  section  4.2)  . 

Split  x  =  (  (Xjj)  i  *  1,  ....  k  j  =  1,  ...,ci  )  into  two  groups  y  and  z. 

Let  y  =  (  (xy)  j  =  1,  ...,ci  )  where  i  takes  on  the  k-m  values 

corresponding  to  the  uncontrolled  chains. 

Let  z  =  (  (xy)  j  =  1,  ...,ci  )  where  i  takes  on  the  m  values 
corresponding  to  the  controlled  chains. 


I 


| 

E* 


(5.2.10) 


Using  this  split  of  x  into  y  and  z,  e(x)x  splits  as  follows. 

e(x)x  =  C|(y,z)y  +  e^y.z)* 

where 

et(yfz)  is  a  k  by  k-ra  matrix 
e^z)  is  an  k  by  m  matrix 

Plugging  equation  5.2.10  into  equation  5.2.8  gives 

yL(y^)cl(y^)y  =  ^(y.zXfCy.z)  -  e2(y,z)z)  (5.2.1 1) 

If  the  k-m  by  k-m  matrix  g-^fy.zje^.z)  is  invertible,  then  from  equation  5.2.11  we  get 

y  =  [^(y.z^^y.z)]  *  gi(y^)(f(y,z)  -  e2(y,z)z)  .  (5.2.12) 

Note:  when  the  k-m  by  k-m  matrix  g^ly.zje^y.z)  is  singular  or  nearly  singular,  the  nonlinear 
inverter  will  produce  very  large  signals  which  will  typically  be  unacceptable. 

If  any  controller  uses  the  m  controls  to  control  the  m  z -chains  to  constant  values,  then  z  =  0 
and  equation  5.2.12  reduces  to 

y  =  [gJ-(y,z)e,(y,z)  ]  1  r^y.ztfly.z)  .  (5.2.13) 


If  f(y,z)  is  polynomial  in  the  entries  of  y,  then  let  y  be  the  column  vector  containing  the  p 
monomials  in  the  entries  of  y  on  which  f  depends.  In  this  case  f  can  be  rewritten  as 

f(y,z)  =  f(z)y  (5.2.14) 

where  f(z)  is  a  k-m  by  p  matrix  of  coefficients.  Equation  5.2.13  now  becomes 

y-  [ff1(y,z)e1(y,z)]  *  [g^y.zjfUjJy  .  (5.2.15) 

Now  consider  the  special  case  where  g^Cy.z)  and  e^y.z)  do  not  depend  on  y. 

Let  R(z)  =  [g^z^z)]"’  [^(z)f(z)]  (5.2.16) 


(5.2.17) 


Using  R,  equation  5.2.15  can  be  written  as 

y  =  R(z)  9 

Equation  5.2.17  is  a  polynomial  ODE  in  y  ,  for  each  fixed  value  of  z  .  The  equilibria  of  these 
ODE  are  given  by  the  solution  of  the  polynomial  equations  : 

0=  [^(zjfo)]*  .  (5.2.18) 

The  entries  of  the  k-m  by  p  matrix  R(z)  are  referred  to  as  the  stability  parameters  since  they 
determine  the  stability  of  the  uncontrolled  chains  of  integrators. 
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5  J  Aircraft  Complementary  Dynamics 


Introduction 


The  6DOF  aircraft  equations  of  motion  can  be  written  in  the  form: 

e(x)x  =  f(x)  +  g(x)h(x,u) 

where  x  =  (V,(3,a,<i>,0,vi\<X>,9) . 

Equation  5.3.1  can  then  be  put  in  the  form 

et(z)y  +  e2(y,z)z  =  f(z)y>  +  g(z)  h(y,z,u) 


(5.3.1) 


(5.3.2) 


where  y  =  (V.'F)  ,  z  =  (|3,a,<i>t0,<I>,©)  ,  and  y  = 


Equation  5.3.2  is  in  the  form  required  for  the  calculation  of  the  complementary  dynamics 
discussed  in  subsection  5.2. 

We  will  derive  the  expressions  for  et(z) ,  c^(y,z) ,  f(z) ,  and  g(z)  .  The  equations  of 
motion  of  an  aircraft  are  simplest  in  body-axis  coordinates  (U,V,W,P,Q,R,<1>,0)  ,  so  we  will 
start  with  them  ,  then  change  coordinates  to  x  =  (V.p.a.O,©,'?,^,©)  .  The  reason  for  chang¬ 
ing  to  these  new  coordinates  is  that  the  complementary  dynamics  are  simpler  then. 


» 

$ 

£ 


& 


% 

l 

! 
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Dimensionless  equations  of  morion  in  (U,V,W,P,Q,R,C>,6)  coordinates 

In  order  to  make  the  equations  of  motion  easier  to  manipulate,  we  will  take  various 
groups  of  expressions  which  appear  together  and  give  them  a  combined  name.  The  most 
natural  way  to  do  this  is  to  form  various  dimensionless  groups  of  expressions. 

(N.B.  :  to  distinguish  dimensionless  parameters  from  physical  ones,  we  could  have  put  tildas 
over  each  dimensionless  variable  in  the  remaining  part  of  the  report.  We  decided  not  to  for 
notational  ease.) 

In  the  following  pages,  it  will  be  shown  how  to  make  the  equations  of  motion  presented 
in  section  3.3  take  on  the  following  dimensionless  form: 

Dimensionless  force  equations  (forces  divided  by  mg)  : 


V  [ul  ,  TP' 

V  =  Q  V  +  1^  +  Q  +  Sxyz  h(x.u)  . 

w  LwJ  IrJ 


Dimensionless  moment  equations  (moments  divided  by  {mg  d^})  : 

?  [p1  Ip' 

LQ  =  nimo  Q  +  V2^  +  Q  +giJim  h(x,u)  . 

R  IBJ  LrJ 


(5.3.3) 


(5.3.4) 


Each  term  in  the  above  equations  is  dimensionless.  The  terms  in  these  two  equations 
will  be  described  in  the  next  page. 

In  equations  5.3.3  and  5.3.4  ,  all  variables  were  made  dimensionless  using  the  following 
constants: 

Let  V  =  a/-2!S- 

nom  V  VipS  ' 

Then  define  t^m  ,  (!,„„,  ,  and  I^  in  terms  of  V,^  ,  m,  and  g  ,  the  acceleration  due  to 
gravity. 

V 

*nom  —  ~  •  ^nom  —  •  an<^  4iom  = 


WWW 


To  make  the  notation  more  compact,  we  will  combine  the  dimensionless  chord  and 

mean  aerodynamic  chord  ,  .  L  ,L  .  , 

,  c  - - ~r -  )  with  the  associated 

dnom 


,  ,  wing  span 
wingspan  (  b  =  — ~  r — 

%om 


aero  functions  as  follows: 


C*yz  = 


C,(a,p) 

Cy(a,{3) 

Cr(a,P) 


*y*fo« 


=  'A 


Cx,(a,P)  C^(a,p)  C^a.P) 
Cy,(a,P)  CyQ(a,p)  Cy/tt^) 
C^(a,p)  C^P)  C^p) 


b  0  O' 

0  c  0 

0  0  b 

^-lmn  ~ 


b  0  0 
0  c  0 
0  0b 


C,(a,p) 

CJa,?) 

C„(a,P) 


b  0  O' 

0  c  0 

0  0  b. 

C,p(cc.P)  C^CouP)  C,,(a,P) 


b  0  0 

0  c  0 

0  0  b. 

Note  that  b  and  c  are  themselves  dimensionless  now,  so  the  above  expressions  are  still  dimen¬ 
sionless. 


All  speeds  are  divided  by  V^,,,  (  eg  U  =  s-EPe<*.  )  and  the  derivatives  are  taken 


with  respect  to  t  = 


time 


tj*om 


P,  Q,  and  R  are  angular  rates  multiplied  by  tnom  .  Ixx  ,  Iyy  , 


etc.  are  moments  of  inertia  divided  by  I^n 


0  R  -Q 

-sin(0) 

^xx  — ^xy  — ^xz 

Q  = 

-R  0  P 

•  1<t>e  - 

cos(0)sin(O) 

>  ^mo 

~^xy  ^yy  ~^yz 

o 

°\ 

O 

,cos(0)cos(O) . 

~^xz  — ^yz  ^zz 

gxyz  = 


C*T  C,, 

Cyr  *-y«  ^y». 


CjT  CL  CL  cl 


.  and  Slmn  = 


b  0  0 
0  c  0 
0  0b 


Ck  CU  Cii 


On*  Qro*t 

Q*x  Qi»  Qi* 


. 


Change  coordinates  from  (U,V,W)  to  (V,p,a)  using  the  following  formula  : 


•1 

u 

cos(a)cos(p) 

V 

=  Vlpa 

where  lpa  = 

sin(p) 

LwJ 

,sin(a)cos(p) . 

Differentiating  gives 


V 

• 

V 

V 

=  Lt  Lq 

p 

w 

a 

where 


►  » 
cos(a)cos(P)  -cos(a)sin(P)  -sin(a) 

10  0 

L,= 

sin(P)  cos(P)  0 

and  Lq  = 

0  V  0 

_sin(a)cos(P)  -sin(a)sin(p)  cos(a) 

0  0  Vcos(P) 

Change  coordinates  from  (P,Q,R,<I>,0)  to  (<i>,0,vF,<I>,©)  using  the  following  formula  : 


v 

Q 

=  ^2 

0 

'F 

where 


L2  = 


1  0  -sin(0) 

0  cos(d>)  cos(0)sin(d>) 
0  — sin(O)  cos(0)cos(d>) 


Differentiating  gives 


where 


P 

<i> 

Q 

=  L 2 

0 

+  Lj 

0 

R 

*  4 

53 


0  0  0  0  0  -cos(©) 

Lj  =  0  — sin(<I>)  cos(9)cos(<l>)  <i+  0  0  -sin(©)sin(d>)  © 

0  -cos(<I>)  -cos(©)sin(<l>) .  0  0  -sin(©)cos(<I>) 

The  expressions  for  P,  Q,  and  R  can  also  be  put  into  &  to  give 

O  =  QjO  +  02©  +  Oj^P 


(S.3.7) 


where 


b  0 

o‘ 

0 

-sin(<!>) 

-cos(<I>) 

o,= 

0  0 

1 

,  02=  sin(<D) 

0 

0 

0  -1 

0. 

[cos(<I>) 

0 

0 

0  cos(©)cos(<I>)  -cos(0)sin(<I>) 

O3  =  -cos(©)cos(<i>)  0  -sin(©) 

cos(0)sin(<I>)  sin(Q)  0 


Substituting  all  this  into  equations  5.3.3  and  5.3.4  gives  the  equations  of  motion  in 
(V,P,a,<X>,0,*P,<I>,©)  coordinates. 


Force  equations  : 

y 

Li  Lo  P 

a 


=  OV  lpa  +  l«e  +  +  VCxyifg(t  L2  9  +  g xyi  h(x,u) 


(5.3.8) 


Moment  equations  : 

*  r  * 

4> 

Ijno  ^2  9  +  L«2  © 

S'  xp 


<t>  <j> 

=  Ol^Lj  $  +  V2^  +  VC^  L2  0  +  glmn  h(x,u)  (5.3.9) 
'F  T 


}•* 

'< 

*;< 

,*!* 


1 


a 

& 
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Expressions  for  et(z) ,  e-^y.z) ,  and  f(z) 


Equations  5.3.8  and  5.3.9  can  be  rewritten  as 


\UU  0*,  1  j  .  .  .  { 

0,c  i^j  i  =  «p,a.®.e.®.6)  +  g(i)h(y,z,u) 


(5.3.10) 


To  fit  f(z)  into  the  space  on  the  page,  split  it  into  its  first  three  columns  and  its  last  three 
columns. 


Let  f(z)  =  j?l23(z)  f456(z)] 


In  the  expression  for  f^z) ,  use  the  following  notation 


&  <j>  o 

Q4  as  +  fljQ  ,  L5  =  L2  9  ,L^  —  L2  9  ,  Lj  *  Lj  0 


d>  <D 

+  ,  Lj  0  =  L5+l*e'F  .and  L*  0  =  L6 +  L/F  .(5.3.10b) 

xp  xp 


Substituting  equation  5.3.10b  into  equations  5.3.8  and  5.3.9  gives 

Q^pa+Qyz**1^  O3 

f|J3<2)  =  C^L,  n<i™!«9+nJin,,L5-UL1 


(5.3.11) 


^xyz  ^3^Pa"^-xyz^Q«^<t>0 


f456(z)  Ci,,,,, 


(5.3.12) 


,i?  *iVlV 


r«wr< 


wvfi1 


.r 


In  equation  5.3.10  ,  g(z)  is  given  by 


g(z)  = 


g*yi(M 

glmn(P»®) 


and  an  example  of  what  h(y,z,u)  can  be  is  given  by 

Thrust 


h(y,z,u)  = 


mg 

V2  sin(8,-  ~) 

V2  sin(8e-a) 
V2  sin(8r-p) 


Equation  5.3.10  can  be  rewritten  as 


where 


°3 

°3 


VL3  Ojjj 
°3x2  W-4 


=  f(2) 


1 

V 

* 

V2 

VT 


+  g(z)h(y,z,u)  (5.3.13) 


^3  = 


L'*'  J 

-cos(a)sin(p)  -sin(a)cos(P) 

1  0 

cos(P)  0 

and  L4  = 

0  cos(<X>) 

-sin(a)cos(P)  cos(a)cos(P)  _ 

0  -sin(d>) 

so 


e,(z)  = 


l|ia  °3 

O3  W**© 


and  e^Cy.z)  = 


vl3  °3x2 
03*2  ^mo^4 


This  is  now  in  the  form  of  equation  5.3.2  so  we  can  apply  the  results  of  subsection  5.2  . 
First  multiply  equation  5.3.2  by  g-*-(z) 

(where  g-**(z)  is  a  2  by  6  matrix  which  satisfies  g^z)  g(z)  =  0  ,  see  section  4.3  ). 

g-L(z)el(z)y  +  g1(z)e2(z)z  =  ^(z^y  (5.3.14) 


y  =  [Wzje^z)]  gHzXffz)?  -  e2(y,z)z)  . 


(5.3.15) 


Equation  5.3.15  holds  in  general.  When  z  =  0  ,  equation  5.3.15  reduces  to 

y=  |gL(z)ei(z)j  g1(z)f(z)y  . 


(5.3.16) 


and  the  second  and  third  columns  of  f(z)  go  to  zero,  leaving  only  columns  1,  4,  5,  and  6  so 
equation  5.3.16  reduces  to 


- 

1 

V 

=  jg^(z)e|(z)]  Vcz)  £ft(z)  f4(z)  f5(z)  f6(z)  ] 

V2 

VY 

(5.3.17) 


Equation  5.3.17  represents  the  complementary  dynamics.  It  is  a  system  of  2  quadratic 
O.D.E.’s  with  the  following  coefficients  (where  subscripts  represent  the  powers  of  V  and 
respectively) : 


aoo  aM  au  a^ 
boo  b20  bn  bo2 


[gJ“(z)el(z) J  V(z)  [f^z)  f4(z)  fj(z)  f6(z)  ]  .  (5.3. IS) 


The  inverse  of  gx(z)el(z)  needed  in  equation  5.3.18  can  be  computed  as  follows.  Start  by 
partitioning  the  2  by  6  g^(z)  matrix  into  four  parts. 


g^z)  = 


•  ■> 

St  Sb 

Sc  gd. 


(5.3.19) 


where  g,  ,  gb  .  gc  .  and  g<j  are  each  size  1  by  3  .  Using  this  notation. 


gJ-(z)e1(z)  = 


Si'!*pa  Sb^mo  l<t>© 
gcJpo  gdl  mo  ^<t>0 


(5.3.20) 


Since  gJ-(z)e1(z)  is  a  2  by  2  matrix,  it  is  trivial  to  invert. 


For  aircraft  whose  controls  have  left/nght  symmetry,  the  elevator  and  thrust  only  affect 
the  first,  third,  and  fifth  rows  of  (5.3.1)  while  the  rudder  and  aileron  only  affect  the  second, 
forth,  and  sixth  rows  of  (5.3.1)  .  The  result  is  that  the  g(x)  matrix  then  has  the  following 
simplified  form: 


g(z)  = 


CxT 

0 

0 

0 

cy. 

0 

°y\ 

0 

0 

0 

bCu 

a 

0 

bCu 

r 

0 

CQn*i 

0 

0 

0 

(5.3.21) 


In  this  simpler  situation,  the  transpose  of  (4.5.8)  becomes  : 


transpose 


(^(z)]  = 


Cx/V  C* 
0 


l  CyiCn*.  Cy*.Cn*r 

b  Cu  C*  -  Cu 


at  r  a 


l  cxTcz*t 

c  ~  <*,  Cnv,. 


(5.3.22) 


1  cy»,ci«r  ^i^u, 

b  <W  Cj, 

or,  in  case  the  denominators  in  5.3.22  are  zero,  use 


.vAMUlAAA  X. 


b(CliCn4-CuC„4) 


transpose 


[s^-w]  - 


•  r  t  • 


^yi  ~  ^y»  Q14 

*r  •  •  r 


cy».C|*,  Cy»,Cl«. 


c^C^-C^) 


c(^'XrQn«> 


CXtpLt~  CXtQz4 


In  the  further  simplifying  situation  where  the  thrust  is  along  the  x  axis  (QT  =  0  =  C^) 
and  the  ailerons  only  produce  a  rolling  moment  ( Cyt  =  0  =  ),  Equation  5.3.22  reduces  to  : 


g^z)  = 


0  10  0  0 


0  0  10- 


(5.3.23) 


Using  this  gx(z)  along  with  a  symmetric  aircraft. 


( I*y  =  0  =  Iyz  ,  Cy(a,0)  =  C,(a,0)  =  Qa.O)  =  0  ) 


with  no  dynamic  aero  coefficients. 


and  flying  with  no  sideslip  and  wings  level 


(  P  =  0  =  <D) 


rj  v.  w-u  vv  ) 


rjirjw'jvyy 


reduces  equation  5.3.18  to  the  following  simple  form: 


*00  *20  au  *02  r  ,  i-i  ,  r  _  .  _  i 

.boo  b a,  b„  b^J  =  |ffL«ei<*>J  ^  [  fi<z>  W  f5<z>  f6(z)J=  (5124) 


Cz-C, 


cos(9) 

sin(a) 


m  p 

Vnu 

• 

sin(a) 


^  [cos(Q)  0  sinCQ)]!^!^ 
cC,^  sin(a) 


cos(y)_i 

_ 

cosCQiljj+sinC©)!^ 


where  y  =  ©  -  a  when  <J>  =  0  =  P  .  The  eight  parameters  in  the  above  system  typically 
have  the  following  signs. 


+  -  0  + 

0  0-0 

which  results  in  a  globally  stable  phase  portrait  in  the  physical  (  V  >  0  )  half  of  the  (  V  ,  4* ) 
phase  space. 


5.4  Canonical  Forms  for  the  Complementary  Dynamic  Equations 


We  have  seen  that  the  complementary  dynamic  equations  are  of  the  form: 


• 

V 

aoo  +  a20^2  +  all^/  + 

y 

boo  +  b20V2  +  bnW  +  b02'F2 

(5.4.1) 


These  equations  describe  the  dynamical  behavior  of  the  aircraft  when  the  control  inputs 
are  used  to  keep  the  direction  of  the  velocity  vector  (a, (5)  and  the  attitude  (<X>,0)  constant 
Stable  equilibrium  solutions  represent  steady-state  flight  conditions  -  we  are  therefore 
interested  in  determining  when  equations  5.4.1  admit  stable  equilibria.  Furthermore,  we  would 
like  to  characterize  the  dynamical  behavior  of  the  solutions  away  from  equilibria  to  describe 
the  aircraft’s  transient  behavior.  The  transient  behavior  of  the  model  should  reflect  some  of 
the  flying  qualities  of  the  aircraft  during  precision  pointing  or  tracking  tasks.  We  will  show  in 
the  next  section  how  to  compute  parameters  that  quantify  the  transient  behavior  in  terms  suit¬ 
able  for  flying  quality  evaluation.  Specifically,  we  will  construct  Liapunov  functions  that 
describe  the  stability  of  the  equilibria  and  provide  a  bound  on  the  settling  time  to  equilibrium 
from  a  non-equilibrium  condition. 

Before  constructing  the  Liapunov  functions,  we  will  transform  the  physical  (V.'F)  coordi¬ 
nates  into  a  new  set  of  coordinates  (x,y)  that  describe  the  system  more  economically.  The  new 
coordinates  will  be  related  to  the  original  (V.'F)  coordinates  by  an  invertible  linear  transfor¬ 
mation.  In  the  new  coordinates,  the  Liapunov  functions  have  an  especially  simple  form  that 
renders  the  flow  in  the  phase  space  easy  to  visualize. 


The  coordinate  transformation  is  developed  through  a  computational  algorithm  based  on 
manipulation  of  quadratic  forms.  Note  that  equations  (5.4.1)  can  be  written: 


V=  [l  V  t] 


aoo 


0 


n  &u 

0  a20  ~ 

n  311 
0  ~2~  302 


,  4>=  [l  V  H*) 


boo  0  0 

bjj 

0  *>20  ~y~ 

bj  i 

o  -r  boz 


1 

V 

I'i'J 


‘  (5.4.2) 


In  this  way  the  complementary  dynamics  are  defined  in  terms  of  the  quadratic  forms  A 
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and  B  represented  by  the  3x3  matrices  above  whose  entries  are  simple  functions  of  the 
coefficients  ay  and  by.  For  convenience  of  notation,  we  write  the  above  equations: 


V  =  A(V.'F)  ,  4'  =  B(V,4'). 


(5.4.2*) 


We  need  to  know  how  the  equations  above  transform  when  the  coordinates  of  the  phase 
space  undergo  a  linear  transformation.  Choose  coordinates  (x,y)  related  to  the  (V.'F)  coordi¬ 
nates  by  an  invertible  transformation  F: 

f"  fi  H  =  [V1  .  (5.4.3) 

fn  %  y.  v 


'l'  ►  ^ 

V=  [l  X  y]x  x  .  [l  x  y]§  x 

lyJ  ly. 


(5.4.4) 


where 


1  0  *  1  0 


A=  0  F  A0  F  ’  B=  OF  B  0  F 


1  0  1  0 


(5.4.5) 


In  the  convenient  notation,  V  =  A(x,y) ,  V  =  B(x,y).  To  eliminate  V  and  T,  use  G  = 


fxl  fell  8l2  V 

lyJ  =  fe.2  822]  1*  • 


(5.4.6) 


to  find 


? 

I 


(5.4.7) 
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The  formula  in  equation  5.4.7  describes  the  complementary  dynamics  in  the  new  set  of 
coordinates  (x,y).  Any  property  of  the  system  5.4.7  that  remains  valid  under  linear  transforma¬ 
tions  is  also  a  property  of  the  original  system  5.4.1,  and  vice-versa.  In  particular,  the  stability 
of  either  of  these  two  systems  implies  the  stability  of  the  other. 

The  reason  for  introducing  new  coordinates  is  to  represent  the  dynamic  equations  using 
quadratic  forms  that  are  as  simple  as  possible.  We  think  of  the  change  of  variables  as  a 
transformation  (A,B)  -->  (A,B)  used  to  reduce  the  number  of  independent  variables  required 
to  represent  the  quadratic  forms.  There  is  another  transformation  useful  for  the  same  purpose; 
we  discuss  it  next. 


Suppose  a  system  is  represented,  as  in  5.4.7,  by  equations  of  the  form: 


=  H 


A(x,y) 

B(x,y) 


(5.4.8) 


for  some  constant  2x2  matrix  H.  Pick  a  nonsingular  2x2  matrix  K,  and  define 
A  =  kuA  +  k12B  ,  B  =  k2iA  +  k^B.  The  same  equation  can  then  be  written 


A(x,y) 

B(x,y) 


(5.4.9) 


Transformations  of  this  type,  which  do  not  involve  a  change  of  coordinates,  may  also  be  used. 

We  use  both  types  of  transformation  in  deriving  canonical  forms  for  the  complementary 
dynamic  equations. 


Theorem  5.4-1:  Suppose  the  complementary  dynamic  equations  5.4.1  admit  some  real  equili¬ 
brium  points  and  are  nondegenerate  (nondegenerate  means  that  at  least  one  of  the  two 
matrices  A  and  B  is  rank  three,  and  that  any  nontrivial  linear  combination  of  them  is  rank  at 
least  two  -  this  is  a  generic  condition).  Then  there  is  a  set  of  coordinates  (x,y),  related  to  (V, 
'P)  by  a  linear  transformation,  for  which  the  complementary  dynamic  equations  are  in  one  of 
the  two  forms: 
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(5.4.10) 


6)  =  h 


1  -  x2  -  y2 
-2xy 


(5.4.11) 


Proof:  Start  with  A  and  B  represented  as  in  5.4.2,  assume  that  A  is  rank  three  (if  A  is  rank 
two,  interchange  A  and  B).  Consider  the  matrix  S  defined  by 


0  0  0 

S  =  booA  —  aooB  —  0  Su  Sj2  • 

0  Sj2  S22 


(5.4.12) 


The  matrix  S  is  symmetric  and,  because  the  equations  are  nonsingular,  rank  2.  We  have 
assumed  that  some  real  equilibria  do  exist,  so  the  nonzero  eigenvalues  of  S  have  opposite 
signs.  Therefore,  for  any  real  c,  there  is  a  3x3  matrix  Fc 


1  0  0 
0  fcll  fc21 
0  fc!2  fc22 


(5.4.13) 


$ 

I 


g 

I 
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such  that 


0  0  0 
FctSFc=  0  0c 
0  c  0 


(5.4.14) 


Pall  the  matrix  on  the  righthand  side  of  5.4.J4  ,Bj.  .Because  of  the  assumed  nondegen¬ 
eracy  of  the  equations,  Fc  can  be  chosen  so  that 


a'oo  0  0 


A'  =  fJaF,  =  I  0  a'. 


(5.4.14’) 


11 

T~  »02 


has  the  property  that  none  of  the  diagonal  entries  are  0. 


From  the  matrix  A’  defined  in  5.4.14’  subtract  — — Bi  to  obtain 

2c  1 


nn  0  0 


0  a  90  0 


(5.4.15) 


0  0  a'n 


Now  by  scaling  x  and  y  independently,  and  by  proper  choice  of  c,  A(  and  B,  can  be 
transformed  simultaneously  to 


1  0  0  0  0  0 

A  =  a<x)  0  e  0  ,  B  =  0  0  -l 
0  0  f  0-10 


where  e  and  f  are  either  1  or  -1.  Because  we  have  assumed  real  equilibria,  e=-l  or  f=-l.  The 
expression  is  symmetric  in  x  and  y,  so  only  two  distinct  cases  arise.  These  are  realized  by 
choosing  e=-l,  and  taking  f=l  or  f=-l. 
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In  the  first  case,  where  f=I,  the  dynamic  equations  are  reduced  to  the  form  of  5.4.10. 
This  is  the  case  where  there  are  two  real  equilibrium  points:  those  corresponding  to  (1,0)  and 
(•1,0)  in  the  (x,y)  coordinates.  In  the  second  case,  where  f=-l,  the  dynamic  equations  are 
reduced  to  the  form  of  5.4. 11.  This  is  the  case  where  there  are  four  real  equilibria:  those 
corresponding  to  (1,0),(0,1),(-1,0),  and  (0,-1)  in  the  (x,y)  coordinates.  The  proof  is  complete. 


In  the  original  (V.'F)  coordinates  the  complementary  dynamic  equations  were 
parametrized  by  8  independent  variables  ay ,  by.  By  the  transformations  shown  above,  the 
number  of  independent  parameters  in  the  equations  can  be  reduced  to  the  4  entries  of  a  con¬ 
stant  matrix  H.  The  H  matrix,  and  the  matrix  transformation  from  the  (V.'f')  coordinates  to  the 
(x,y)  coordinates,  are  all  that  is  required  for  a  complete  analysis  of  the  complementary 
dynamic  equations.  Expressions  for  these  matrices  can  be  computed  from  the  construction 
given  in  the  proof  of  Theorem  5.4-1. 


In  the  next  section  we  construct  Liapunov  functions  for  the  two  canonical  forms  of 
Theorem  5.4-1. 


£ 


Si  The  Liapunov  Functions 


In  the  last  section  we  showed  how  the  complementary  dynamic  equations  could  be 
transformed  to  new  coordinates  in  which  fewer  parameters  are  required  to  specify  the  system. 
One  advantage  of  these  new  coordinates  is  that  they  make  it  easy  to  write  down  Liapunov 
functions  for  the  stable  equilibria  (when  stable  equilibria  exist).  In  this  section  we  show  how 
these  Liapunov  functions  are  constructed  and  explain  what  they  can  tell  us  about  flying  quali¬ 
ties.  In  one  of  the  cases  discussed  below,  a  Liapunov  function  is  used  to  show  global  stability 
of  the  model  in  the  V,*?  phase  plane.  Some  of  the  quantities  associated  with  the  construction 
could  be  used  to  measure  flying  quality. 


We  have  seen  that  the  complementary  dynamic  equations  in  the  coordinates  (x.y)  look 


like 


X 

htI  hl2 

1  +  ey2  -  x2 

y. 

" 

h21  h^ 

-2xy 

(5.5.1) 


where  the  components  hlllhl2,h2i,h22  of  the  matrix  H  are,  for  each  set  of  values  a,p,0,<t>,  a 
set  of  constants  determined  by  the  nonlinear  aircraft  model.  In  the  case  where  there  are  four 
real  equilibria  in  the  V.'F  plane  the  parameter  e  takes  the  value  -1,  while  in  the  case  of  two 
real  equilibria  e  takes  the  value  +1.  This  representation  only  applies  when  there  do  exist  real 
equilibria.  We  do  not  consider  the  case  where  there  are  no  real  equilibria. 


The  first  step  in  constructing  the  Liapunov  functions  is  to  take  the  inner-product  of  the 
above  equation  with  the  vector  (1  +  ey2  -  x2  ,  -2xy)  to  find 


(1  +  ey2  -  x2)x  -  2xyy  =  [(1  +  ey2  -  x2)  ,  -2xy] 


hu  h12 

1  +  ey2  -  x2 

h2,  h22 

-2xy 

(5.5.2) 


The  resulting  expression  on  the  righthand  side  of  the  equation  will  have  a  definite  sign  (except 

at  equilibria)  whenever  the  symmetric  matrix  S  =  -j(H  +  HT)  is  positive  or  negative  definite. 

Let  us  suppose  for  the  moment  that  S  is  positive  definite.  The  lefthand  side  of  equation  5.5.2 
is  in  each  case  an  expression  easy  to  represent  by  the  derivative  of  a  simple  function.  The 
level  sets  of  such  a  function  are  transverse  to  the  flow  and  can  be  used  to  characterize  the  sta¬ 
bility  of  the  equilibria. 
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Of  course,  there  are  matrices  H  for  which  the  quadratic  form  above  is  indefinite.  When  S 
is  indefinite,  it  may  be  difficult  to  determine  much  of  value  from  the  functions  constructed 
below.  Still,  our  experience  in  working  with  the  F-14  model  suggests  that,  for  a  wide  variety 
of  flight  conditions,  the  H  matrix  is  definite,  so  that  the  analysis  below  is  relevant. 

CASE  1:  e=+l  .  . 


In  this  case  there  are  two  real  equilibria:  one  at  (1,0)  and  the  other  at  (-1,0).  For  this 
configuration  there  are  three  possibilities  for  the  stability  of  the  equilibria: 

case  1  -  (1,0)  stable  and  (-1,0)  unstable 
case  2  -  (1,0)  unstable  and  (-1,0)  stable 
case  3  -  (1,0)  and  (-1,0)  saddles. 


Our  construction  will  provide  a  Liapunov  function  in  the  first  two  cases  where  one  of  the 
equilibria  is  stable. 


Define  the  complex  number  z  =  x  +  iy,  and  think  of  the  transformed  V.'F  plane  as  the 
complex  z  plane.  The  function 

|2 

(5.5.3) 


F(z)  =  i-^ — £■ 


1 0  +  z)  | 


is  the  one  we  want  to  analyze.  Except  for  the  values  0,1, «> the  level  sets  of  this  function  are 
circles  in  the  plane  who^e  centers  lie  on  the  real  axis.  The  value  <»  is  realized  only  at  the 
point  z  *  -1  and  0  is  realized  only  at  z  =  +1.  The  1 -level  set  is  the  imaginary  axis  x  =  0. 

From  the  geometry  of  the  level  sets  of  the  function  F,  it  is  clear  that  the  point  z  =  + 1  is  a 
global  attractor  for  the  entire  z-plane  whenever  it  can  be  shown  that  the  function  F  restricted 
to  the  solution  curves  of  the  system  is  decreasing.  Then  F  is  a  Liapunov  function  for  the  flow. 
There  is  an  easy  condition  on  S  to  check  that  will  tell  us  if  F  is  decreasing: 

Lemma  5.5-1:  If  the  matrix  S  is  positive  definite,  then  the  function  F  is  a  Liapunov  function 
for  the  system  5.5.1. 


Proof:  In  terms  of  x  and  y. 
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F(x  +  iy) 


(1  -  x)2  +  y2 
(1  +  x)2  +  y2 


(5.5.4) 


and  so 

-£-F(x  +  iy)  =  - - — 3"2-«l  +  y2  -  x2)x  -  2xyy) .  (5.5.5) 

at  ((1  +  x)2  +  y2)2 

At  every  point  other  than  z  =  -1  the  sign  of  this  expression  is  opposite  to  that  of  the  quantities 
in  equation  5.5.2.  But  if  S  is  positive  definite,  the  sign  of  the  quantities  in  5.5.2  is  positive 
(except  at  z  =  +1  where  the  expressions  vanish).  It  follows  that  at  every  point  other  than  an 
equilibrium  value,  the  trajectories  of  the  solutions  to  the  differential  equations  are  such  that  F 
is  decreasing  along  them.  Then  (1,0)  is  a  globally  stable  equilibrium,  while  (-1,0)  is  a  source. 
The  proof  is  complete. 

Observe  that  the  open  half-plane  x>0  coincides  with  the  set  of  values  (x,y)  for  which 
IFkl.  In  the  case  where  H  is  positive  definite,  we  can  compute  a  bound  on  the  rate  at  which 
the  trajectories  starting  in  this  half-plane  converge  to  the  equilibrium. 

Lemma  5.5-2:  Suppose  S  is  positive-definite,  and  that  F0,  the  value  of  F  at  a  state  (x0,y0)  at 
time  t=0,  is  smaller  than  1.  Let  X^fS)  denote  the  smaller  eigenvalue  of  S.  Then  for 
(x(t),y(t))  along  the  solution  trajectory: 

F(t)  <  F 0e'4X^(S)t  (5.5.6) 


Proof:  By  the  hypothesis  on  S 


_d_p  <  ~4\nin(S)(l  +  X2  +  y2)2 
dt  <  ((1  +  x)2  +  y2)2 

<  -^(SjF 


(5.5.7) 


because,  for  all  (x,y). 
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(5.5.8) 


r  .  «  -  »>2  +  y 

(1  +  x)2  +  y2 


<  (i  +  x2  +  y2)2 

((1  +  X)2  +  y2)2 


It  follows  that  F  is  approaching  0,  its  value  at  the  equilibrium  point  (1,0),  at  least  as  fast  as 
the  solution  to  the  bounding  equation 

=  -^WS)F  (5.5.9) 

But  the  solution  for  this  equation  in  F  satisfying  the  initial  condition  F(0)  =  F0  is  exactly  the 
righthand  side  of  5.5.6.  The  proof  is  complete. 


In  the  case  where  S  is  negative  definite,  the  same  function  F  is  used  to  show  that  (-1,0) 
is  globally  stable,  while  (1,0)  is  the  source.  In  either  case,  the  eigenvalue  of  S  having  the 
smaller  magnitude  gives  some  quantitative  estimate  of  how  stable  the  stable  equilibrium  is.  If 
iXnunl  is  close  to  zero,  for  example,  the  system  will  be  slower  to  converge  than  if  it  is  large. 
Also,  it  is  more  likely  when  1  I  is  small  that  a  small  perturbation  to  one  or  more  of  the 
aircraft  model  parameters  could  produce  instability  in  a  nominally-stable  case.  In  general, 
when  that  smallest  magnitude  eigenvalue  of  S  is  close  to  0,  the  pilot  will  have  to  wait  for  a 
while  before  the  aircraft  settles  into  a  steady-state  condition  when  he  commands  a  constant 
turn  at  a  fixed  0-<I>  attitude.  This  parameter  should  also  affect  the  quality  of  the  ride  during 
more  dynamic  maneuvers. 

When  the  unique  physical  (V  >  0)  equilibrium  is  stable,  the  analysis  of  this  case  has  a 
simple  interpretation.  To  reach  a  desired  trim  condition,  the  pilot  need  only  get  the  angles 
a,|3,e,<I>  to  the  correct  values  and  wait.  If  his  speed  is  too  slow,  his  speed  will  increase  to  the 
proper  value.  If  too  fast,  it  will  slow  down.  Likewise,  the  heading  rate  4*  will  find  its  equili¬ 
brium  value  and  stay  there.  The  stability  of  the  unique  equilibrium  insures  that,  if  the  angles 
are  right,  the  vehicle  will  naturally  take  care  of  the  rest. 

CASE  2:  e»-l 


In  this  case  there  are  four  real  equilibria:  (1,0), (0,1), (-1,0), (0,-1).  By  a  general  result  of 
Kukles  and  Casanova  (see  (KCj  or  theorem  7  of  [Cop])  two  of  these  are  saddles,  one  is  a 

70 


***  i.+  i 


sink,  and  the  last  one  is  a  source.  The  sink  and  source  must  be  opposite  each  other  (i.e.  nega¬ 
tives  of  each  other)  so  there  are  four  distinct  possibilities:  any  one  of  these  four  points  could 
be  the  unique  attractor.  We  consider  only  the  case  where  (1,0)  is  the  attractor,  so  that  (-1,0)  is 
the  source  and  the  points  (0,1)  and  (0,-1)  are  saddles.  The  other  three  cases  can  be  converted 
into  this  one  by  a  suitable  linear  transformation,  so  there  is  no  need  to  analyze  the  other  cases 
separately. 


The  function  we  consider  for  tnis  situation  is  the  polynomial 


F(x,y)  =  x(l  -  y2  -  yx2) . 


(5.5.10) 


Lemma  5.5-2:  If  the  matrix  S  is  positive  definite,  then  the  function  F  is  a  Liapunov  function 
for  the  system  5.5.1. 


Proof:  It  is  easy  to  see  that 


-j-F(x,y)  =  (1  -  x2  -  y2)x  -  2xyy  . 
dt 


(5.5.11) 


Comparing  the  righthand  side  of  this  expression  with  the  quantities  in  equation  5.5.2,  we  find 
that  it  is  positive  (except  at  the  equilibria)  whenever  S  is  positive  definite.  So  F  is  an  increas¬ 
ing  function  along  the  solution  curves  for  5.5.1.  The  proof  is  complete. 


To  see  what  this  means,  consider  the  0-level  set  for  the  function  F.  The  0-level  set  is  the 


union  of  the  y-axis  (x=0)  and  the  ellipse  C  defined  by  1  -  y2  -  yx2  =  0.  The  two  saddles 


(0,1)  and  (0,-1)  lie  in  this  set  (they  are  in  fact  the  points  of  intersection  of  the  line  and  the 
ellipse),  while  the  two  other  equilibria  are  contained  in  the  region  bounded  by  the  ellipse. 
The  point  (1,0)  is  a  local  maximum  for  the  function  F,  and  the  point  (-1,0)  is  a  local 


minimum. 


If  the  matrix  S  is  positive  definite,  then  the  function  F  increases  along  the  solution  curves 
of  the  differential  equation.  Therefore,  any  point  inside  the  region  bounded  by  the  y-axis  and 
the  right  half  of  the  ellipse  C  will  be  moved  along  a  trajectory  that  stays  inside  that  region 
and  approaches  the  point  (1,0).  The  region  bounded  by  the  y-axis  and  the  right  half  of  C  is  a 
domain  of  attraction  for  the  point  (1,0).  The  mirror  image  region  in  the  left  half-plane  is  a 
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region  of  repulsion  for  unstable  point  (-1,0). 

The  saddles  at  (0,1)  and  (0,-1)  make  the  global  properties  of  the  flow  unstable.  If  the  air¬ 
craft  state  gets  into  the  unstable  flow  region  of  one  of  these  saddles,  however,  the  pilot  will 
have  to  change  his  aircraft’s  a,P,0,<I>  values  or  fall  into  an  unstable  spin  or  speed  condition. 
The  a,p.0,<!>  values  that  give  rise  to  these  saddles  represent  potentially  hazardous  conditions, 
while  the  region  bounded  by  the  0-level  set  in  the  right  half-plane  defines  the  region  where  it 
is  safe  to  fly. 

There  are  two  physical  equilibria  in  this  case,  at  most  one  of  them  is  stable.  Even  if  there  is 
a  stable  one,  the  simple  strategy  of  setting  the  angles  a,p,0,<I>  to  the  proper  values  and  wait¬ 
ing  will  not  always  work.  If  the  pilot  has  just  completed  a  maneuver  that  has  left  his  V  and  4* 
states  in  a  bad  spot  (too  near  the  unstable  equilibrium)  and  then  attempts  to  keep  the  angles 
fixed  at  the  correct  values  for  the  stable  equilibrium,  he  could  find  himself  in  a  divergent 
speed  condition  or  an  unstable  spin.  If  his  initial  V  and  4*  are  in  a  good  spot  (near  enough  to 
the  stable  equilibrium),  on  the  other  hand,  he  will  be  fine.  Trying  to  end  maneuvers  at  equili¬ 
bria  like  these  could  be  a  hazardous  undertaking. 


SECTION  6:  NONLINEAR  FLYING  QUALITIES 


In  previous  sections  we  discussed  nonlinear  models  of  aircraft,  the  trim  set,  and  the  tech¬ 
nique  of  dynamic  inversion.  These  three  topics  are  basic  to  our  understanding  of  nonlinear 
flying  qualities,  which  we  discuss  in  this  secdon. 

We  begin  with  a  discussion  of  two  idealized  types  of  parameters:  commanded-dynamic 
parameters  and  complementary-dynamic  parameters.  Both  sets  of  functions  are  computed 
directly  from  the  nonlinear  aircraft  models  and  they  quantify  important  physical  properties  of 
the  vehicle’s  behavior  during  flight.  The  commanded-dynamic  parameters  measure  the  maneu¬ 
verability  and  controllability  of  the  aircraft  in  the  three  angle-rate  degrees  of  freedom  (P,Q,R). 
The  commanded-dynamic  parameters  are  discussed  in  subsection  6.1  below. 

Subsection  6.2  covers  the  complementary-dynamic  parameters.  Complementary-dynamic 
parameters  can  be  used  to  construct  Liapunov  functions  for  stability  analysis  of  dynamic 
inversion  controllers.  The  best  results  we  have  so  far  are  for  the  a,fi,0,<I>  inversion,  for  which 
we  have  derived  explicit  time  and  space  bounds  for  the  nonlinear  dynamic  trajectory  of  the 
vehicle  moving  towards  a  trim  condition. 

Besides  these  two  types  of  idealized  parameters,  we  have  found  several  criteria  for  super- 
maneuverable  vehicles  flying  along  trajectories  where  aerodynamic  forces  and  inertial  terms 
simultaneously  play  an  important  role.  We  have  two  main  results  here 

1)  a  simple  aerodynamic  criterion  for  smoothness  of  the  aerodynamic  loading  during  rapid  a 
variation 


*>] 


ft 


2)  control  design  criteria  for  coordinated  flight  during  highly-dynamic,  simultaneous  roll- 
pitch-yaw  maneuvers. 


These  two  results  are  discussed  in  subsections  6.3  and  6.4 

The  last  two  topics  of  this  section  arc  two  ideas  that  were  under  development  at  the  end 
of  our  program.  The  first  is  a  flying-quality  metric  for  nonlinear  aircraft  models  that  could  be 

73 


evaluated  using  the  coordinated-flight  U,P,Q,R  dynamic- inversion  controller  discussed  in  sec¬ 
tion  6.4.  The  second  concerns  the  stability  and  controllability  of  the  aircraft  during  maneuvers 
involving  extreme  angular  rates.  The  parameters  defined  in  this  section  quantify  the  effect  of 
the  dynamic  derivatives  on  the  rotational  energy  and  angular  momentum  of  the  vehicle. 

Many  of  the  ideas  below  were  inspired  by  analysis  of  maneuvers  like  those  described  in 
section  7,  which  follows  this  one.  The  maneuvers  we  have  looked  at  were  generated  by  the 
batch  version  of  our  flight  simulation  program  using  various  dynamic  inversion  strategies. 
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6.1  Commanded  Dynamic  Parameters 


Commanded  dynamic  parameters  quantify  the  pilot’s  direct  command  authority  over  the 
state  of  the  aircraft  by  use  of  controls.  Examples  of  parameters  from  the  MIL-F-8785C 
specifications  of  this  sort  are  (by  section  number): 

3.2.2. 1  Short-period  response 

3.22.2  Control  feel  and  stability  in  maneuvering  flight  at  constant  speed 

3.2.3. 3.1  Longitudinal  control  in  catapult  takeoff 

3.2.3.4  Longitudinal  control  in  landing 

3. 3.2.6  Turn  coordination 

3.3.4.(all)  Roll  control  effectiveness 

3.3.5  Directional  control  characteristics 

3.3.7  Lateral-directional  control  in  crosswinds 

3.3.8  Lateral-directional  control  in  dives 

3.3.9  Lateral-directional  control  with  asymmetric  thrust 


3.4.2. 1.3  Stall  prevention  and  recovery 

All  these  criteria  concern  command  response  —  the  response  of  the  vehicle  in  degrees  of  free¬ 
dom  that  the  pilot  is  intentionally  changing  through  direct  control.  In  contrast,  there  are  cri¬ 
teria  concerning  the  response  of  the  aircraft  in  degrees  of  freedom  that  the  pilot  is  not  trying 
to  change  during  the  course  of  a  maneuver  (e.g.  3. 2. 1.1. 2  Pitch  control  force  variations  during 
rapid  speed  changes). 
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In  an  effort  to  identify  useful  new  flying  quality  parameters,  we  defined  some  functions 
of  the  aircraft  state  and  the  nonlinear  aircraft  model  that  quantify  fundamental  limitations  on 
dynamic-maneuver  command  responses.  These  functions  quantify  the  command  authority 
available  to  the  pilot  to  control  vehicle  dynamics  during  maneuvers.  Below  are  some  exam¬ 
ples  of  these  functions  applied  to  the  (P,Q,R)  command-response.  Dimensionless  coordinates 
are  used  throughout  this  section  (see  section  5.3  for  the  conventions). 

Pitch-rate  Control 


The  primary  pitch-rate  effector  in  a  basic  aircraft  is  the  elevator.  In  dimensionless  coordinates 
(assuming  symmetric  aircraft,  no  dynamic  derivatives,  and  a  simplified  elevator  model): 

IyyQ  =  (R2  -  P2)!,*  +  RP(In  -  I„)  +  V2^  +  V^sintfe  -  a)  (6.1.!) 


The  derivative  of  Q  has  an  inertial  term,  an  aerodynamic  term,  and  a  control  term  depending 
on  the  elevator.  A  parameter  which  measures  the  basic  Q  command  effectiveness  is 


V2^ 

*yy 


(6.1.2) 


This  parameter  must  be  sufficiently  large  at  low  speed  to  assure  adequate  control  of  the  pitch 
axis  for  take-off,  landing,  and  near-stall  maneuvers. 

Another  potentially  useful  function  is  the  dynamic  pitch-control  ratio,  defined  as: 


VQc  = 


(R2  -  P2)I„  +  RP(IU  -  I„)  +  V2^ 


V^, 


m* 


(6.1.3) 


This  function  measures  the  ratio  of  the  uncommanded  portion  of  Q  to  the  magnitude  of  the 
elevator  authority.  In  a  region  of  the  state  space  where  this  ratio  is  too  large,  the  pilot  will 
have  trouble  maintaining  attitude  control.  For  example,  if  the  ratio  is  bigger  than  1,  the  pilot 
cannot  even  control  the  sign  of  Q.  Some  degradation  of  pitch  control  during  high-a  or  rolling 
maneuvers  is  expected,  quantifies  the  extent  of  degradation. 

If  R  and  P  are  0,  reduces  to  the  ratio  of  Cm  to  Cm< .  In  this  case,  the  parameter  is  a 
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function  of  a  alone  (assuming  coordinated  flight)  --  it  could  be  used  to  determine  if  the  vehi¬ 
cle  has  adequate  elevator  authority  for  a  desired  pullup  maneuver. 

Roll-rate  and  Yaw-rate  Control 

For  a  symmetric  aircraft,  the  P  and  R  degrees  of  freedom  are  most  naturally  analyzed 
together.  The  basic  nonlinear  equations  look  like: 


I** 

-I«  I** 


0  R  -Q 
Q  -P  0 


0  -lxz 

0  lyy  0 

-i«  o  Ia 


p 

C|' 

cu. 

c'..' 

Q 

R 

+V2 

4- 

s 

s 

V2sin(5r-  P) 


In  coordinated  flight  (P  =  0),  symmetry  causes  the  aerodynamic  functions  Cj  and  Q  to 
vanish,  resulting  in  a  simplification  of  the  righthand  side.  When  the  angular  velocity  vector 
(P,Q,R)  is  small,  roll/yaw  command  authority  is  approximated  by  the  matrix: 


^PRar 


^xx  “^Xl 

-i 

cT 

rh z 

Cn*. 

(6.1.5) 


The  size  of  the  minimum  singular  value  c^n^pR*-)  is  a  bound  on  the  angular-rate  authority 
available  in  some  direction  in  this  two-degree-of-freedom  subspace.  The  singular  vectors 
(input  and  output)  associated  with  this  singular  value  should  be  considered  as  well.  The  func- 
tion  cfmin(TlpR„)  depends  on  the  speed,  a  and  P,  and  the  vehicle  model.  The  larger  it  is,  the 
better  the  pilot  can  control  P  and  R  independendy.  It  can  be  thought  of  as  a  lateral  - 
directional  version  of  the  parameter  defined  earlier. 

For  more  dynamic  but  still  coordinated  maneuvers,  the  ratio  is  defined  to  be  the 
vector: 
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(6.1.6) 


The  two  entries  in  this  vector  represent  the  size  of  the  aileron  and  rudder  settings  required  to 
match  the  inertial  1  and  n  components  of  the  torque.  At  states  where  these  values  are  large, 
the  pilot  must  use  large  rudder  and  aileron  commands  to  maintain  or  reduce  the  sizes  of  P  and 
R.  Where  these  values  are  too  large,  the  spin  condition  might  be  beyond  the  pilot’s  control. 


The  parameters  rij^  and  Vp^  can  be  defined  for  uncoordinated  flight  as  well.  All  that 
is  required  is  to  include  the  aerodynamic  terms  in  Vp£ff. 


6.2  Complementary  Dynamic  Parameters 


We  define  complementary  dynamics  to  be  those  dynamics  of  the  vehicle  associated  with 
states  which  the  pilot  is  not  controlling  directly.  Depending  on  the  pilot’s  task,  the  states 
associated  with  the  complementary  dynamics  may  vary.  We  consider  the  following  MIL-F- 
8785C  specifications  (listed  by  their  MIL-F-8785C  section  numbers)  to  be  associated  with 
complementary  dynamic  phenomena: 


3.2.1. (all) 

Longitudinal  static  stability 

3.2.2. 1.3 

Residual  oscillations 

3.2.2.2 

Control  feel  and  stability  in  maneuvering  flight  at  constant  speed 

3.3.1.(all) 

Lateral-directional  mode  characteristics 

3.3.2. 1 

Lateral-directional  response  to  atmospheric  disturbances 

3.3.5. 1 

Directional  control  with  speed  change 

3.3.9.3 

Transient  effects 

3.4.2. 1.2 

Stall  characteristics 

3.4.2.2 

Post-stall  gyrations  and  spins 

3.4.2.2.1 

Departure  from  controlled  flight 

3.4.3 

Cross-axis  coupling  in  roll  maneuvers 

3.4.4 

Control  harmony 

3.4.5 

Buffet 

3.4.11 

Direct  force  controls 

3.5.5. 1 

Failure  transients 

3.6.3 

Transients  and  trim  changes 

All  these  specifications  concern  the  response  of  the  aircraft  to  effects  other  than  those  directly 
commanded  by  the  pilot.  Included  are  transient  response  (in  degrees  of  freedom  other  than 
those  commanded),  disturbance  and  failure  response,  dynamic  cross-coupling,  and  stability 
and  damping  of  dynamical  modes. 

We  have  identified  some  parameters  based  on  the  nonlinear  models  that  describe  one 
type  of  complementary  dynamic  behavior,  we  call  them  complementary  dynamic  parameters. 
Our  approach  makes  use  of  the  nonlinear  inversion  method  discussed  in  section  5. 


In  the  basic  example  where  the  pilot  has  four  control  inputs  to  work  with,  he  will  be  able 
to  command  independently  only  four  degrees  of  freedom,  leaving  two  degrees  of  freedom 
constrained  by  their  relation  to  the  others.  The  motion  of  the  aircraft  in  those  two  remaining 
degrees  of  freedom  is  determined  by  the  equations  of  motion  and  the  pilot's  command  inputs, 
but  the  control  the  pilot  has  in  those  two  dimensions  is  indirect  -  a  side-effect  of  the  direct 
control  exerted  in  the  other  four  dimensions.  The  dynamic  behavior  of  the  aircraft  in  the  four 
dimensions  where  the  pilot  has  control  is  (within  the  limits  of  control  effectiveness  and 
neglecting  disturbances)  determined  by  the  control  inputs,  so  the  pilot  can  directly  influence 
dynamics  there.  Once  the  pilot’s  choice  is  made,  the  dynamics  in  the  other  two  dimensions 
are  completely  determined  by  the  aerodynamic  and  physical  properties  of  the  aircraft.  These 
residual  or  indirectly  controlled  dynamics  will  vary  from  one  aircraft  to  another  in  a  way  that 
can  be  computed  from  the  nonlinear  models.  From  these  dynamics,  flying  quality  parameters 
can  be  computed. 

To  compute  these  parameters,  we  make  assumptions  about  the  strategy  chosen  by  the 
pilot  to  control  the  aircraft.  For  illustration,  suppose  the  pilot  uses  his  four  controls  to  keep 
the  direction  of  his  velocity  vector  (a  and  P)  and  the  direction  of  the  gravity  vector  (0  and  <j>) 
fixed.  This  strategy  might  be  used  by  the  pilot  to  execute  a  steady  turn.  As  was  shown  in 
section  5.3,  the  complementary  dynamics  for  V  and  'F  are  given  by  the  equations: 

V  =  ago  +  a^V2  +  auV'F  +  agj'F2  (6.2.1) 

*  =  bgg  +  ^V2  +  bnVY  +  boz'F2 

The  stability  of  these  complementary  dynamic  equations  should  be  highly  correlated  with 
the  ride  quality  during  the  steady  turn.  If  the  system  is  very  stable,  the  ride  should  feel  very 
steady  (good  flying  qualities  for  tracking  purposes,  for  example),  but  if  it  is  only  marginally 
stable  the  ride  quality  may  feel  unsteady  or  even  oscillatory.  In  the  worst  case,  if  the  system 
is  unstable,  the  aircraft  might  experience  a  divergent  speed  condition  (compare  with  MIL-F- 
8785C  section  3.2. 1.1)  or  an  unstable  spin  condition  (compare  with  MIL-F-8785C  sections 
3.4.1  and  3.4.2).  Analysis  of  the  stability  of  this  system  was  performed  in  subsection  5.5  of 
the  previous  section.  We  derived  Liapunov  functions  to  describe  the  transient  dynamics  for 
V  and  4/,  and  we  computed  a  bound  on  the  settling  time  from  a  non-equilibrium  condition. 
In  one  of  the  two  cases  analyzed  in  that  subsection  (corresponding,  perhaps,  to  a  very  extreme 
turn  or  a  turn  at  high  a),  we  found  that  the  stability  of  the  tum  might  depend  on  the  values  of 


V  and  4*  as  the  turn  is  started.  The  functions  and  parameters  of  that  section,  and  other  param¬ 
eters  computed  in  a  similar  fashion,  seem  very  important  in  assessing  how  well  aircraft  enter 
and  execute  turns. 

By  fixing  other  combinations  of  states  (or  some  four-dimensional  subspace  of  the  state 
space),  other  sets  of  analytic  equations  can  be  obtained.  From  these  equations,  other  comple¬ 
mentary  dynamic  parameters  can  be  computed. 


6J  :  Lift  to  Drag  Ratios  of  Aircraft  with  Smooth  Lift  Curves 


Many  books  on  aerodynamics  give  estimates  for  the  lift  and  drag  of  a  wing  (or  aircraft  with 
zero  control  deflection).  These  estimates  are  typically  of  the  form 


CL  =  fi(AR,a) 


Q> -  cDo  * 


(C  Q2 
f2(AR) 


where  AR  is  the  aspect  ratio  of  the  wing  (or  aircraft)  . 

Taking  the  ratio  of  the  equations  6.3.1  and  6.3.2  gives 

CD  ~  Cd»  f,(AR,a) 

CL  =  f2(AR) 

If  the  lift  were  linear  in  a  ,  then  we  would  get 

Cl  =  fi(AR.ot)  =  Cl^(a-<x0) 


(6.3.1) 


(6.3.2) 


(6.3.3) 


(6.3.4) 


and  equation  6.3.3  would  be  linear  in  a 

fj(AR.a)  C^o-Oq) 

f2(AR)  f2(AR) 


(6.3.5) 


For  lift  curves  that  are  not  linear  in  a  ,  we  would  still  expect  equation  6.3.3  to  be  fairly  linear 
in  a  for  small  a  . 

For  several  aircraft  with  smooth  lift  curves  (no  stall  discontinuities)  such  as  the  F-4,  F-14  and 
F-15,  we  have  examined  the  low  speed  lift  and  drag  curves  measured  during  wind  tunnel  tests 
and  found  that  the  expression  on  the  left-hand  side  of  equation  6.3.3  was  very  nearly  deter¬ 
mined  by  the  following  simple  and  interesting  form  : 


(6.3.6) 


''Do  7C 

- — - =  tan(a  -  clq)  (  0  ^  a  £  —  radians  ) 

Cl  2 

This  formula  appears  in  "USAF  Datcom  Methods  Handbook  for  Double  Delta  Wings”  . 
The  point  of  giving  it  hear  is  to  show  how  well  it  fits  the  F-4,  F-14,  and  F-15  data  for  angles 
of  attack  ranging  from  0  to  beyond  90  degrees. 

In  (6.3.6),  0^  is  defined  to  be  the  angle  of  attack  at  which  the  drag  is  minimum.  Note  that 
tan(a  -  a^)  is  nearly  linear  in  a  -  Oq  ,  over  a  fairly  large  range  of  a.  -  clq  . 

For  the  F-4  and  F-15,  AR  is  approximately  2.9  ,  while  on  the  F-14  (which  had  wing  sweep 
set  at  22°  )  AR  was  approximately  7.3;  yet  equation  6.3.6  still  holds.  This  indicates  that 
even  though  the  numerator  and  denominator  in  equation  6.3.3  may  each  depend  on  AR  ,  their 
ratio  does  not  depend  on  AR  for  these  aircraft. 


An  interesting  interpretation  of  equation  6.3.6  can  be  made  when  a  change  of  coordinates  is 
made  from  wind  axes  to  body  axes. 


Let  ft  =  a  -  otQ  ,  then 


cos(a) 

sin(a) 


■ 

-sin(a) 

cos(a) 


(6.3.7) 


so 

cx  =  -cos(a)CD  +  sin(a)CL  =  -cos(a)(CD  -  CdP  +  sin(a)CL  -  cos(a)CDo  (6.3.8) 


Plugging  6.3.6  into  6.3.8  gives 


Cx  =  -Co,  cos(ft) 


(6.3.9) 


If  equation  6.3.6  were  exact  then 

I  Cx  I  £  (  0  S  a  £  radians  ) 


(6.3.10) 


For  the  F-4,  F-14,  and  F-15  data  shown  in  figures  6.1,  6.2,  and  6.3  at  the  end  of  this  subsec¬ 
tion,  Cd,  is  around  .02  and 

ICxl£3Ci^  (0£a£y  radians)  (6.3.11) 

Since  is  so  small,  the  error  in  the  CD  and  CL  data  from  which  Cx  was  calculated 
may  account  for  much  of  the  discrepancy  between  (6.3.10)  and  (6.3.11)  . 

Equation  6.3.11  shows  us  that  at  low  speed  there  is  much  less  aerodynamic  force  in  the  x 
direction  than  in  the  z  direction  (  V2  C*  is  small  at  any  angle  of  attack  ).  Consequently,  the 
only  significant  force  the  pilot  feels  in  the  x  direction  is  due  to  the  throttle.  This  is  true  no 
matter  how  rapidly  a  is  varying.  Therefore,  it  would  seem  that  aircraft  with  small  magnitude 
have  better  flying  qualities  during  rapid  a  maneuvers  than  those  with  Cx  of  large  variation. 
Note  also  from  the  plots  that  Cx  decreases  (it  is  negative)  almost  monotonically  with  a  ,  for 

0  £  a  £  “  radians  .  This  may  be  of  some  use  in  using  it*  at  the  percussion  point  to  "meas- 

ure"  a  (using  a  lookup  table  for  Cz  ). 

The  F -15  data  (at  0  sideslip)  came  from  the  tests  conducted  on  a  13-percent  scale  model  of 
the  F-15  S/MTD  configuration  in  the  NASA-Langley  30x60  ft  low-speed  wind  tunnel.  Test 
No.  489  ,  run  216  ,  no  canard  ,  configuration  85.000  ,  between  10  July  and  27  July  1985. 

The  F-4  data  (at  0  sideslip)  and  the  F-14  data  (at  22°  wing  sweep,  sideslip  *  0°  and  20°) 
were  taken  from  [MMTJ]. 

This  report  cited  reference  [Ang]  for  wind  tunnel  data  for  the  F-4  and  several  references  con¬ 
taining  selected  data  for  the  F-14  from  wind  tunnel  tests  conducted  during  15  March  -  16 
April,  1971,  in  the  NASA  Ames  Research  Center  12  ft  pressure  tunnel  and  during  August, 
1971,  in  the  NASA  Langley  Research  Center  30x60  ft  (full  scale)  tunnel. 
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6.4  Highly  Dynamic  Phenomena 

For  a  conventionally  configured  aircraft  at  medium  to  high  speed,  the  largest  forces  that 
can  act  are  aerodynamic.  Even  the  most  powerful  aircraft  have  thrust-to-weight  ratios  only 
slightly  greater  than  1,  while  maneuvers  reaching  levels  close  to  10  g  for  short  periods  of  time 
are  not  uncommon.  For  coordinated  flight,  these  high-g  maneuvers  are  achieved  by  entering 
high-angle-of-attack  regions  where  the  aerodynamic  lift  and  drag  forces  are  strongest.  One 
example  of  this  type  of  maneuver  is  the  diving-tum  maneuver  shown  later  in  section  7.3. 

The  largest  force  during  this  diving-tum  maneuver  is  the  8-g  normal  acceleration  encoun¬ 
tered  during  the  start  of  the  dive.  In  less  than  3  seconds,  alpha  rises  from  close  to  zero  to 
nearly  1.1  radians.  The  vehicle  decelerates  rapidly,  despite  the  fact  that  it  is  diving  at  full 
throttle,  because  the  drag  term  is  so  large.  The  elevator  is  saturated  at  1  radian  deflection  (our 
assumed  saturation  value)  to  maximize  the  pitch-rate  during  this  period  -  the  high  pitch-rate 
begins  while  the  roll-rate  is  still  large  from  the  initial  banking  phase.  The  large  pitch-rate 
increases  alpha  rapidly  while  the  airspeed  is  still  near  500  ft/second,  large  lift  and  larger  drag 
forces  are  the  result. 

An  important  feature  of  this  high-g  diving  maneuver  is  the  dominance  of  the  quadratic 
inertial  terms  in  the  nonlinear  state  equations.  That  is,  the  angular  rates  P,  Q,  and  R  become 
so  large  in  magnitude  that  the  state  derivatives  are  heavily  influenced  by  the  product  terms 
PQ,  QR,  PW,  etc.  A  primary  effect  of  these  large  angular  rates  can  be  seen  in  the  basic 
lateral  velocity  equation: 


V  =  -RU  +  PW  +  cos(0)sin(4>)  +  V2Cy(a,(3)  +  direct  rudder  acceleration  (6.4.1) 


This  simple  model  (no  dynamic  derivatives,  only  rudder  direct  forces)  illustrates  the  point.  If  a 
maneuver  involves  turn-coordination  at  high  roll  rate,  the  magnitudes  of  the  terms  -RU  and 
PW  become  comparable  with  the  gravity  term,  and  much  greater  than  the  direct  rudder  force. 
In  our  sample  maneuver,  the  roll  rate  P  reaches  a  value  greater  than  two  radians/second  and 
the  pitch  rate  exceeds  1  radian/second  in  the  early  pan  of  the  maneuver.  The  maneuver  begins 
at  at  a  speed  of  roughly  500  ft/second  (size  2  in  the  dimensionless  units  of  equation  6.4. 1 )  and 
at  an  altitude  of  about  5000  ft  (near  sea  level,  but  at  least  2700  ft).  The  inertial  acceleration 
Ny  experienced  by  the  pilot  is  the  negative  of  the  sum  [V2Cy(a,P)  +  direct  rudder  force],  and 
this  is  a  small  term  because  the  turn  is  coordinated.  The  gravity  term  never  gets  any  larger 
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than  1,  while  each  of  the  terms  -RU  and  PW  get  bigger  than  1. 

To  keep  P  small,  the  yaw-rate  R  must  vary  in  such  a  way  as  to  keep  V  small.  The  con¬ 
troller  used  to  generate  the  diving  turn  allowed  arbitrary  commanded  P,  and  chose  R  to 
satisfy: 


R  =  [PW  +  cos(0)sin(<J>)  +  V2Cy(a,P)  +  direct  rudder  acceleration]/U  (6.4.2) 

The  control  strategy  based  on  the  dynamic  inversion  for  U,P,Q,R  with  R  chosen  to  kee  3 
small  is  called  coordinated-flight  U,P,Q,R  inversion.  This  controller  was  used  to  execute  the 
diving-turn  maneuver  of  section  7.3,  it  seems  a  very  promising  approach  for  supermaneuver- 
able  vehicle  control.  If  the  lateral-directional  control  authority  is  adequate  to  generate  the 
necessary  P  and  R  values  independendy,  it  should  work  for  very  extreme  combined  roll- 
pitch-yaw  maneuvers. 

For  this  control  approach  to  work  in  practice,  the  allowed  bandwidth  of  the  roll-rate  com¬ 
mands  must  be  limited  to  a  region  where  the  rudder  is  able  to  generate  enough  R  to  track  the 
right-hand  side  of  equation  6.4.2.  As  a  part  of  this  restriction,  consideration  must  be  give,  to 

w 

a  because  P  is  multiplied  by  —  =  tan(a).  The  factor  tan(a)  implies  that  at  large  angles  of 
G)p(a) 

attack  the  ratio - of  the  bandwidths  must  be  smaller  than  at  low  a  for  coordinated  rolls. 

toR(a) 

Equation  6.4.2  can  be  used  to  quantify  that  requirement:  for  acceptable  high-alpha  roll-rate 
response,  the  bandwidths  coR(a)  and  0)p(a)  should  be  related  by  an  inequality  of  the  type: 

coR(a)  >  Ka>p(a)tan(a)  (6.4.3) 

for  some  constant  K  which  is  (presumably)  larger  than  1.  This  inequality  can  be  viewed  as  a 
design  requirement  for  acceptable  flying  qualities  at  high  a  (better  flying  qualities  being  asso¬ 
ciated  with  larger  values  of  K),  or  it  can  be  taken  as  a  rule  for  determining  the  maximum 
allowable  bandwidth  of  P  as  a  function  of  a  and  the  bandwidth  of  R.  Alternatively,  it  can  be 
viewed  as  a  requirement  on  rudder  size  needed  to  provide  coordinated  flight  at  high  a.  The 
parameter  Vp^  defined  in  section  6.1  must  also  be  considered  to  determine  the  achievable 
performance  for  this  control  approach. 
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Along  the  same  lines,  it  is  possible  to  derive  criteria  for  dynamic  control  of  a  in  terms  of 
Q  control  During  coordinated  flight,  P  and  R  do  not  contribute  directly  to  a,  but  Q  does. 
Assuming  V  »  0: 


a  =  Q  +  [Ucos(0)cos(4>)  +  Wsin(0)  -  UN,  +  WNJ/V2  (6.4.4) 


Consider  a  maneuver  that  requires  keeping  a  large  for  an  extended  period  of  time.  From 
6.1.1,  it  is  clear  that  if  P  and  R  are  kept  small,  the  sign  of  Q  is  determined  by  the  sign  of 
at  values  of  a  where  1 I  >  I  CL  I  .  If  is  negative  with  a  large  magnitude  in  the  desired 

Tl 

a  range,  then  Q  will  be  decreasing  rapidly  so  long  as  a  is  kept  large.  Once  Q  becomes 
sufficiently  small,  the  sign  of  a  in  equation  6.4.4  becomes  negative  and  the  desired  large  a 
condition  is  lost.  From  this  analysis,  we  can  see  one  of  the  advantages  of  a  thrust-vectoring 
capability.  When  the  factor  C,,,  dominates  CL,  for  values  of  a  in  a  desired  range,  thrust  vec- 

toting  can  maintain  Q  so  long  as  V2  is  small  enough.  In  such  a  situation,  it  is  possible  to 
maintain  a  large  a  condition  for  extended  periods  of  time  without  generating  P  and  R  (i.e., 
turning)  to  keep  Q  in  a  compatible  range. 


6.5  Dynamic  Flying  Quality  Metrics 


One  approach  to  evaluating  flying  qualities,  given  a  nonlinear  aircraft  model,  is  to  use  a 
metric  on  the  space  of  maneuvers.  This  idea  has  been  proposed  by  others,  we  have  chosen  not 
to  develop  such  an  approach  during  this  contract.  A  metric  approach  could  be  carried  out 
based  on  ideas  developed  during  this  contract;  one  way  to  do  this  is  described  briefly  below. 

Begin  by  choosing  a  (small)  number  of  equilibrium  conditions  for  the  aircraft  model. 
Dynamic  maneuvers  begin  at  various  equilibrium  conditions;  the  equilibrium  points  chosen 
should  be  representative  of  flight  conditions  where  dynamic  maneuvers  usually  start.  We 
could  describe  the  set  of  points  chosen  by  their  equilibrium  values  of  a,(i,0,  and  O  (these 
coordinates  are  the  easiest  to  use  with  our  equilibrium  computing  program).  For  each  point 
we  then  compute  the  fixed  values  of  V.'F,  and  the  actuator  inputs.  To  reduce  the  number  of 
cases  considered,  we  could  start  with  coordinated  conditions  where  (5  =  0. 

Next,  for  each  flight  condition,  pick  a  set  of  dynamic  command  profiles  for  U,  P,  and  Q 
as  a  function  of  time.  Many  different  sets  of  profiles  can  be  specified  here  *  the  idea  is  to 
define  a  representative  sampling  of  dynamic  responses  for  as  many  different  maneuvers  as 
possible.  The  aircraft  model  will  be  evaluated  with  respect  to  this  set  of  candidate  maneuvers, 
so  it  would  be  a  good  idea  to  include  profiles  that  are  characteristic  of  maneuvers  actually 
used  during  flight. 

Each  U,P,Q  profile  chosen  will  then  be  used  as  input  to  the  U,P,Q,R  dynamic  inverter. 
The  R  command  needed  by  the  controller  will  be  generated  internally  to  provide  coordinated 
flight  while  the  simulation  tries  to  fly  the  trajectory  defined  by  the  specified  U,P,Q  functions. 
For  each  maneuver,  the  model  will  be  evaluated  according  to  a  collection  of  norms  chosen  to 
reflect  the  quality  of  the  response. 

For  example,  for  each  command  profile,  we  would  find  the  average  error  in  the  com¬ 
manded  U,  P,  and  Q  responses  along  the  trajectory;  as  well  as  the  maximum  error  for  each  of 
these  variables.  Also  important  is  the  determination  of  the  extent  of  control  saturation  (both 
position  and  rate  limits)  during  the  flight  These  measures  of  U,  P,  and  Q  errors  for  each 
maneuver  would  form  a  criterion  for  evaluation  of  the  attitude-rate  response,  a  quantity 
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important  in  dynamic  tracking  or  pursuit  tasks.  By  taking  the  average  value  of  these  errors 
over  a  representative  sample  of  maneuvers,  the  evaluator  finds  an  overall  sense  of  the  quality 
of  response. 

The  physical  factors  about  the  vehicle  reflected  by  the  attitude-rate  response  errors  are 
things  like  the  mass  and  moment  properties,  span  and  chord,  aerodynamic  moment  coefficients 
C|,Cin,Cn,  and  the  actuator  position  and  rate  limitations.  The  aerodynamic  force  terms 
Cj.Cy.Cj  will  probably  not  be  very  important  for  this  criterion.  It  might  work  out  better  to  fix 
the  throttle  setting  at  a  predetermined  value  for  the  entire  maneuver  so  that  power  cycling 
uncharacteristic  of  actual  maneuvers  will  not  occur.  In  that  case,  only  the  profiles  for  P  and  Q 
need  to  be  specified  for  each  run.  It  is  probably  a  good  idea  to  classify  the  responses  accord¬ 
ing  to  whether  or  not  some  saturation  occurred  during  the  simulation,  and  to  quantify  the  ten¬ 
dency  to  saturate  in  an  average  sense.  Some  evaluation  should  also  be  made  of  the  severity  of 
the  effects  of  saturation  (i.e.  whether  instability  results). 

Another  factor  to  consider  is  the  attitude  response  during  the  maneuver,  as  well  as  the 
final  equilibrium  position  and  time  to  settle  at  the  end  of  the  maneuver.  For  simplicity,  we 
might  suppose  that  the  commanded  values  of  P  and  Q  go  to  0  at  the  end.  Some  (presumably) 
short  time  thereafter  the  vehicle  should  show  some  tendency  to  settle  at  a  new  equilibrium, 
which  can  be  considered  the  starting  point  of  a  subsequent  maneuver  not  yet  determined.  The 
length  of  the  path  through  the  attitude  space,  and  the  location  of  and  time  to  settle  at  the  final 
equilibrium  point  are  all  quantities  that  can  be  measured  and  averaged  over  the  maneuver 
space.  The  time  to  settle  at  the  end  of  a  maneuver  is  a  quantity  similar  in  nature  to  our 
complementary-dynamic  parameters  that  we  developed  during  this  program. 

Finally,  we  would  look  at  the  inertial  accelerations  during  the  maneuver  to  assess  the 
environmental  stress  experienced  by  the  pilot  during  the  maneuver.  The  values  of  Nx,  Ny,  Nz 
would  be  computed  along  each  trajectory  and  compared  with  acceptable  ranges.  In  the  basic 
nonlinear  model,  the  formulas  for  these  inertial  accelerations  felt  by  the  pilot  are: 

Nx  =  -(V^^a.P)  +  throttle  acceleration)  (6.5.1) 

Ny  =  -(V2Cy(a,P)  +  rudder  acceleration) 

Nz  =  -(V2Cz(a,P)  +  elevator  acceleration) 
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It  would  also  be  worthwhile  to  compute  the  acceleration  time  derivatives,  or  jerk  components, 
during  each  maneuver.  The  jerk  components  quantify  the  smoothness  of  the  ride  in  a  way  that 
allows  comparison  among  various  models.  The  jerk  components  can  be  computed  using  impli¬ 
cit  differentiation  along  the  trajectory:  for  example,  the  x -component  of  the  jerk  Jx  is: 


-ovvc,(o,p) + v^d + 4£-p>  + 

dt  da  dp  dt 


)  (6.5.2) 


with  similar  expressions  for  Jy  and  Jz.  Some  overall  average  of  the  expression  Jx  +  Jy  + 
for  the  different  command  profiles  might  provide  a  criterion  for  smoothness  of  ride  during 
dynamic  maneuvers.  Note  that  the  jerk  components  depend  explicitly  on  the  partial  derivatives 
of  the  aerodynamic  force  functions  -  these  derivatives  could  be  estimated  numerically  from 
interpolated  table  data. 


6.6  Stability  and  Controllability  of  Rotational  Energy,  Angular  Momentum,  Angular  Rate 


One  idea  that  arose  near  the  end  of  the  program  was  to  identify  quantities  based  on  non¬ 
linear  aircraft  models  that  characterize  the  stability  and  controllability  of  the  rotational  energy, 
angular  momentum,  and  rotation  rate  of  the  vehicle  during  flight.  In  this  section  we  show  how 
such  quantities  can  be  derived  and  show  some  examples.  As  a  corollary  of  this  analysis  we 
derive  formulas,  defined  in  terms  of  aerodynamic  functions  and  the  moment  tensor,  which 
quantify  the  stability  of  the  rotational  dynamics  for  bounded  speed. 

Using  the  nonlinear  equations  that  include  dynamic  derivatives,  we  have: 


V  —  +  ^2Qiyz  +  ^OiyzPQR 
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These  analytic  expressions  for  the  derivatives  of  V,  P,  Q,  and  R  involve  coefficients  that 
depend  in  a  complicated  way  on  the  states  o,|3,0,<|>;  but  their  dependence  on  V,  P,  Q,  and  R  is 
only  quadratic.  There  are  V2  terms  in  the  control  input  (g  . . .  h(x,u))  expressions  associated 
with  the  surface  effectors,  if  we  identified  this  dependence  explicitly  in  the  equations  the  form 
would  remain  quadratic.  From  physical  considerations  the  V2  terms  in  equation  6.6.1  represent 
drag  effects  that  dominate  V  whenever  V  is  large  (relative  to  P,  Q,  and  R).  These  equations 
are  similar  to  the  ones  we  encountered  in  the  analysis  of  the  complementary  dynamic  parame¬ 
ters,  only  in  this  case  there  are  no  assumptions  made  about  any  of  the  states  being  kept  fixed. 
We  expect  they  can  be  analyzed  thoroughly  to  obtain  a  global  stability  result,  but  for  now  we 
do  not  attempt  a  complete  analysis. 

Let  us  suppose  that  the  speed  V  is  bounded,  and  consider  equation  6.6.2.  First,  observe 
that  the  dot  product  of  the  vector  [2P,2Q,2R]  with  both  sides  of  equation  6.6.2  provides  an 
expression  for  the  time  derivative  of  P2  +  Q2  +  R2,  the  square  of  the  rate  of  the  angular 
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speed.  The  right-hand  side  is  then  a  cubic  function  of  the  variables  P,  Q,  and  R  with 
coefficients  depending  on  the  speed  V,  the  states  a,P,0,6,  and  the  control  inputs.  We  have  not 
worked  much  with  this  expression  because  the  other  functions  below  seemed  more  tractable. 
Still,  it  might  be  worth  looking  at  to  determine  the  effectiveness  of  the  control  inputs  to  con¬ 
trol  the  angular  speed  directly. 


Simpler  to  analyze  and  (perhaps)  more  useful  is  control  of  the  rotational  energy  about  the 
center  of  mass.  The  time  derivative  of  the  rotational  energy  is  obtained  by  taking  the  dot  pro¬ 
duct  of  both  sides  of  equation  6.6.2  with  the  vector  [P,Q,R]  after  first  multiplying  by  the  iner¬ 
tia  matrix  The  resulting  expression  on  the  righthand  side  is  now  simply  a  quadratic 
expression  because  [P,Q,R]  £2  =  0.  The  equation  becomes: 


-—[Rotational  Energy]  =  Q  ] 

dt  dt  2  R 


(6.6.3) 
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The  expression  we  want  to  consider  is  the  quadratic  term  [P,Q,R]  Q  that  appears 

[rJ 

on  the  righthand  side.  For  bounded  V,  this  quadratic  term  will  dominate  the  rotational  energy 
derivative  when  the  angular  rate  vector  [P,Q,R]  is  large.  If  the  3x3  matrix 

SimnPQR  =  ~(ClmnpQR  +  CjmPQR)  is  negative  definite  for  all  values  a  and  p,  the  rotational 

rate  must  stay  bounded.  Physically,  this  matrix  represents  an  angular  drag  term  -  its  dominant 
effect  on  the  rotational  energy  at  high  angular  rate  conditions  suggests  it  may  be  highly  corre¬ 
lated  with  flying  qualities  for  extremely  dynamic  maneuvers.  The  sizes  of  the  eigenvalues  and 
the  directions  of  the  corresponding  eigenvectors  in  the  P,Q,R  space  are  reasonable  candidate 
flying  quality  parameters. 


The  total  angular  momentum  squared  is  another  function  of  the  angular  velocity  vector 
that  can  be  analyzed  to  quantify  angular  stability  and  control  authority  during  flight.  Proceed¬ 
ing  as  before,  we  compute 
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~  [Angular  Momentum]2  *  -7-[-7[P,Q,R]I^o 
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(6.6.4) 


In  this  case,  the  cubic  P,Q,R  term  vanishes  because  the  matrix  is  antisymmetric.  Now 

it  is  the  matrix 


AimnPQR  =  y(1moOmnPQR  +  CJ^jpqrIX,)  (6.6.5) 

that  must  be  negative-definite  to  insure  stability.  The  eigenvalues  of  this  matrix  and  the 
corresponding  eigenvectors,  computed  as  a  function  of  a  and  f},  might  also  be  correlated  with 
flying  qualities. 

We  do  not  know  which  function  of  the  angular  velocity,  if  any,  is  the  most  important  to 
a  pilot  during  highly  dynamic  maneuvers.  The  relation  between  the  sizes  of  the  matrices 
simnPQR  811(1  AimnPQR  811(1  the  sizes  of  the  different  torque-generation  control  effectors  could 
have  a  large  impact  on  the  vehicle’s  flyability  in  terms  of  the  angular  rate  controllability. 
Parameters  based  on  all  three  functions  discussed  above  have  a  clear  physical  significance  and 
could  be  relevant  to  flying  qualities.  Also,  there  may  be  a  correspondence  between  these 
parameters  and  the  onset  of  uncontrolled  spin  conditions  that  are  sometimes  associated  with 
high-angle-of-attack  maneuvers. 

We  must  use  caution  when  interpreting  the  results  of  this  subsection.  We  do  not  know 
how  large  the  rotational  rates  must  be  for  the  dynamic  derivative  terms  to  dominate.  The 
maximum  rotational  rates  predicted  by  the  bound  on  rotational  energy  from  the  above  equa¬ 
tions  may  exceed  the  realm  where  our  basic  nonlinear  model  is  valid.  We  have  assumed  a  flat 
earth,  constant  air  density,  and  very  simple  expressions  for  the  aerodynamic  functions.  We  do 
know,  however,  that  nonlinear  models  like  these  are  used  in  practice,  and  that  they  give  rea¬ 
sonable  results  when  applied  to  specific  regions  of  flight.  Within  the  realm  where  these  equa¬ 
tions  accurately  model  the  aircraft  dynamics,  we  expect  the  parameters  above  will  be  worth 
looking  at.  If  we  have  a  better  model,  the  three  angular  velocity  functions  discussed  here  still 
have  a  fundamental  physical  significance  that  might  be  correlated  with  flying  qualities. 
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SECTION  7:  THE  MANEUVERS 

In  this  section  we  present  some  simulated  maneuvers  using  a  subset  of  the  F-14  model 
found  in  [MMTJ].  These  simulations  motivated  some  of  the  mathematical  expressions  that  we 
associated  with  candidate  flying  quality  parameters  in  earlier  sections.  More  importantly,  they 
provide  a  detailed  account  of  the  model  behavior  under  carefully  controlled  conditions.  We 
can,  in  general,  investigate  and  demonstrate  dynamic  flight  characteristics  associated  with  vari¬ 
ous  levels  of  flying  quality  through  simulations  like  these. 


Our  approach  has  been  based  on  analysis  of  the  nonlinear  aircraft  equations,  to  character¬ 
ize  the  dynamic  behavior  of  aircraft  during  flight.  At  the  start,  we  began  with  an  idea 
(dynamic  inversion)  about  how  aircraft  could  be  made  to  fly;  we  tried  out  the  idea  by  using  it 
to  command  simulated  maneuvers.  By  analyzing  the  results  of  the  simulations,  we  were  able 
to  identify  features  of  the  nonlinear  equations  that  could  have  a  significant  influence  on  an 
aircraft’s  dynamic  behavior  during  flight  Some  of  the  things  we  learned  are  mentioned  in  the 
presentation  below,  others  we  have  already  discussed  in  our  earlier  section  on  the  candidate 
parameters. 

Our  analysis  of  the  maneuvers  given  here  is  not  so  complete  as  we  would  like  -  there 
remain  several  questions  about  the  results  that  should  be  investigated  further.  We  will  point 
out  the  unsettled  questions  as  we  go. 


The  three  different  types  of  maneuvers  discussed  here  were  generated  by  our  batch  ver¬ 
sion  of  the  nonlinear  simulation.  The  first  type  is  a  roll  reversal  -  where  the  aircraft  banks  first 
to  the  right,  stabilizes  bank  angle,  then  banks  back  to  the  left.  The  second  type  is  a  barrel 
roll.  The  third  type  is  a  highly-dynamic  diving  tum  chosen  because  of  its  extreme  dynamic 
behavior. 
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7.1  Roll  Reversal 


Three  versions  of  the  roll  reversal  maneuver  are  presented.  All  three  were  generated  by  a 
dynamic  inversion  controller  tracking  open-loop  command  profiles  for  a, (3, 741.  The  structure 
of  this  controller  is  described  in  example  5  of  subsection  S.l.  The  primary  goal  was  to  fly  the 
maneuver  while  keeping  the  flight-path  angle  y  fixed  at  zero.  The  three  different  versions 
involved  three  different  criteria  for  the  normal  acceleration  Nr 

The  easiest  way  to  fly  die  maneuver  is  to  pick  a  desired  p  trajectory  and  then  command  y 
-  0,  P  =  0,  and  a  to  be  what  it  has  to  be  in  order  to  keep  the  speed  nearly  constant  The  plots 
in  Figure  7.1  through  Figure  7.8  show  the  results  for  this  maneuver.  We  call  this  maneuver 
UNLOADEDRHVERSAL.  To  generate  UNLOADEDREVERSAL,  we  chose  a  desired  n  com¬ 
mand  profile: 


For  0  <  t  <  5  seconds  commanded  p.  =  0.0 

For  5  <  t  <  6  seconds  commanded  p.  =  ~  sin2(x(t  -  5.0)/2.0) 

For  6  <  t  <  8  seconds  commanded  p  =  -r 

3 

For  8  <  t  <  10  seconds  commanded  p  =  -j  (1  -  2sin2(x(t  -  8.0)/4.0)) 


For  10  <  t  <  12  seconds  commanded  p  =  - 


n 

3 


The  commanded  values  for  P  and  y  were  fixed  at  0  throughout.  To  derive  the  command 
profile  for  a,  we  used  the  approximation  (good  when  y  is  constant  at  0  and  the  direct  control- 
surface  force  terms  are  small): 


y  =  V2CLcos(p)  -  l 


(7.1.1) 


Note  that  variables  in  the  dimensionless  units  of  section  5.3  are  used  here. 

We  wanted  to  keep  y  constant  at  0,  and  we  chose  to  keep  V  near  a  nominal  dimension¬ 
less  speed  of  2.2  (roughly  560  feet  per  second)  throughout,  so  we  needed  to  find  an  a 
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command  profile  (note  that  0  is  assumed  0)  to  solve: 

CL(a)  =  l/((2.2)2  cos(commanded  p))  (7.1.2) 

For  small  a,  the  lift  coefficient  is  closely  approximated  by  the  relation: 

CL(a)  =  5.9a  +  0.09  (7. 1.3) 

Therefore,  we  defined  the  commanded  a  according  to: 

commanded  a  =  (-0.09  +  l/((2.2)2  cos(commanded  p))/5.9  (7. 1 .4) 

It  was  not  essential  that  we  could  approximate  the  lift  curve  by  a  linear  function  of  a,  any 
invertible  nonlinear  function  would  have  worked  just  as  well.  If  a  wider  range  of  a  were 
involved  in  the  maneuver,  we  would  have  used  a  more  accurate,  nonlinear  approximation. 
Plots  of  the  commanded  a  and  p.  profiles,  with  the  simulated  responses,  can  be  found  in  Fig¬ 
ure  7.3. 

First  consider  the  four  plots  in  Figure  7.1.  Each  plot  shows  the  graphs  of  three  of  the 
states  as  functions  of  time  during  the  maneuver.  The  name  of  the  maneuver  and  the  states 
plotted  appears  at  the  top  of  each  plot.  The  way  to  read  these  chans  having  more  than  one 
state  function  plotted  on  a  single  graph  is  as  follows:  the  first-named  variable  at  the  top  of  the 
graph  r  represented  by  the  continuous  line,  the  second-named  variable  is  represented  by  the 
short-dashed  line,  and  the  third-named  variable  is  represented  by  the  long-dashed  line.  In 
some  cases,  one  of  the  three  variables  plotted  remains  much  smaller  in  magnitude  than  the 
others,  so  that  it  appears  that  only  two  graphs  have  been  drawn.  For  example,  in  the  plot  for 
U,  V,  and  W  in  Figure  7.1  the  graph  of  V  is  so  small  that  it  is  difficult  to  tell  that  it  is  there. 
Some  care  is  required  to  read  these  plots  correctly.  In  some  of  these  plots  there  are  rapid  tran¬ 
sients  during  the  first  2  seconds  -  they  are  associated  with  mismatches  of  initial  conditions  at 
the  start  of  the  simulation  and  should  be  ignored. 

There  are  few  features  worth  noting  in  Figure  7.1.  The  components  U  and  V  remain 
nearly  constant  during  the  maneuver  while  W,  though  small,  increases  by  a  factor  of  2  or  3 
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from  beginning  to  end.  The  behavior  of  W  is  driven  by  the  a  response,  which  can  be  seen  in 
Figure  7.3.  The  P,Q,R  plot  shows  that  a  2-radian/second  peak  in  P  and  very  little  Q  and  R  are 
required  to  produce  the  60  degree  bank  in  one  second.  The  stabilization  into  the  banked-tum 

requires  roughly  ~  second.  The  <t»  response  is  crisp,  with  little  overshoot,  and  9  stays  nearly 
* 

fixed.  4*  moves  from  0  to  roughly  15  degrees  with  the  bank  to  the  right,  then  back  to  0  again 
after  the  reverse'  back  to  the  left.  In  inertial  position  coordinates,  the  path  of  the  vehicle  moves 
to  the  right  about  500  feet  during  the  maneuver  from  t=5  seconds  to  t=12  seconds,  and  is 
about  to  start  heading  back  again  when  the  simulation  ends. 

Figure  7.2  shows  plots  of  the  changes  in  the  inertial  position  coordinates  Y  and  Z,  and  of 
the  Euler-angle  coordinates  0  and  H*  drawn  to  appropriate  scale.  Worth  special  note  is  the 
CHANGE  IN  ALTITUDE  plot,  where  it  is  shown  that  the  altitude  varies  by  les:  than  2  feet 
during  the  entire  time.  Keeping  in  mind  that  the  maneuver  did  not  begin  until  t=5  seconds 
(the  initial  condition  transient  was  a  factor  at  the  start),  we  can  see  that  the  altitude  was  held 
very  nearly  constant  while  the  maneuver  was  performed. 

Figure  7.3  shows  the  comparisons  between  the  open-loop  command  profiles  for  a,(J,Y,|J. 
and  the  values  of  these  functions  during  the  simulation.  After  the  initial  condition  transients, 
all  four  simulated  responses  stay  within  a  milliradian  of  their  commanded  values.  This  was  a 
very  successful  maneuver. 

Figure  7.4  shows  the  simulated  values  of  the  inertial  accelerations,  the  aerodynamic  force 
functions,  a  and  (J,  and  the  speed.  Ny  remains  very  close  to  0,  Nx  changes  by  about  g,  and 
Nz  follows  a  benign  path  between  1  and  2  g.  The  speed  changes  by  less  than  2  percent. 

Figure  7.5  shows  the  control  input  behavior  during  the  maneuver.  Disregarding  the  initial 
transients,  these  profiles  seem  fairly  reasonable.  The  aileron  does  most  of  the  work  in  a  pair 
of  doublets  (peak  aileron  deflection  of  50  degrees)  during  the  two  rolling  periods  (we  have 
assumed  1-radian  ranges  for  each  of  the  surfaces  -  we  would  have  assumed  mote  effective 
surfaces  and  performed  the  same  maneuvers  with  smaller  size  commands  if  we  had  less  sur¬ 
face  position  range.  We  will  deal  with  saturation  phenomena  soon,  but  not  for  the  aileron). 
The  elevator  and  rudder  commands  move  only  slightly,  about  5  degrees.  The  surface  com¬ 
mand  rates  look  reasonable,  there  probably  should  be  a  plot  of  (approximate)  surface  rate 
activity  in  a  more  detailed  version  of  this  simulation.  The  throttle  moves  up  to  a  setting  of 
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Figure  7.5  Control  Activity  During  UNLOADEDROLLREVERSAL 


0.5  g  during  the  initial  loading,  then  moves  back  to  a  slower  setting  when  the  banked  turn  is 
established.  During  the  reversal  it  moves  down  to  0.0,  then  back  to  0.4  g  during  the  loading  in 
the  opposite  direction.  The  variation  of  a  during  the  maneuver  changes  the  drag  (see  Figure 
7.4),  so  some  thrust  activity  is  required  if  speed  is  to  be  kept  nearly  constant  Also,  the  thrust 
command  was  varied  by  the  controller  to  generate  small  corrections  to  the  speed  necessary 
because  of  the  direct  force  contributions  of  the  surfaces,  which  we  neglected  when  we  gen¬ 
erated  the  command  profile  for  cl  The  controller  also  had  to  compensate  for  the  fact  that  the 
a-command  profile  was  generated  by  an  approximation  to  the  lift  curve.  The  magnitude  and 
rate  of  the  throttle  activity  during  this  maneuver  might  be  considered  too  large  -  if  so,  a 
different  command  profile  for  clP.y.H  could  be  chosen  or  a  different  inverter  used  for  control¬ 
ling  the  aircraft  The  throttle  was  simply  doing  what  it  had  to  do  to  allow  tracking  of  the 
a,p,y,n  command  profiles  in  Figure  7.3. 

Figures  7.6  and  7.7  depict  the  eight  parameters  a^by  associated  with  the  complementary 
dynamics.  From  them,  we  computed  that  at  each  time  during  the  maneuver  there  was  a  unique 
physical  equilibrium  value  in  the  complementary  V.^F  phase-plane.  Two  of  the  plots  in  Figure 
7.8  show  the  two  eigenvalues  at  that  equilibrium  (both  are  always  real  for  this  maneuver  - 
note  that  the  equilibrium  is  always  stable),  and  the  other  two  plots  compare  the  simulated 
values  of  and  V2  with  the  equilibrium  values.  The  eigenvalues  could  be  used,  as  discussed 
in  section  5.5,  to  compute  a  settling  time  for  the  V.'F  functions  to  their  equilibrium  values 
once  the  states  a,(),0,<I>  have  stabilized.  We  have  not  performed  this  computation  yet,  but  it 
would  not  be  difficult  (we  have  not  coded-up  the  algorithm  on  the  computer  yet).  During  the 
transients  moves  away  from  the  equilibrium  values  by  a  significant  factor  (about  100  per¬ 
cent  error),  but  then  moves  back  to  equilibrium  very  quickly  when  rolling  stops.  V2,  on  the 
other  hand,  tracks  the  equilibrium  value  very  closely  (after  the  initial  transient). 

For  this  first  maneuver,  the  a,P,Y,H  dynamic  inverter  looks  remarkably  good. 

The  second  maneuver  differed  from  the  first  only  in  the  goals  for  the  normal  acceleration 
during  the  reversal.  We  called  it  LOADEDREVERSAL.  The  objective  was  to  keep  Nz  £  2 
after  the  first  60  degree  bank.  Again,  we  chose  to  use  the  a.p.y.H  controller.  The  command 
profiles  for  3  and  p.  remained  the  same  as  in  the  unloaded  case,  we  changed  the  command 
profiles  for  a  and  y.  Our  new  a  command  profile  was  the  same  as  the  old  one  until  t=8 
seconds,  when  the  reversal  began.  To  keep  Nt  above  2  g  during  the  reversal,  we  chose  to 
keep  the  a  command  constant  at  (-0.09  +  2/((2.2)2  ))/5.9  (compare  with  Figure  7.4)  after  t=8 


Figure  7.8  Eigenvalues  and  ,  V2  Comparisons  During  UNLOADEDROLLREVERSAL 


until  the  end  of  the  simulation.  The  y  command  was  not  easy  to  choose:  we  would  have  liked 
to  keep  y  =  0  throughout,  but  it  is  not  physically  possible  to  keep  both  p  and  y  at  0  while 
requiring  that  Nz  be  at  least  2  g  during  the  reversal.  So,  we  decided  to  ask  for  a  small  amount 
of  y  increase  throughout  the  reversal,  the  plot  for  the  command  profile  is  shown  in  Figure 

7.11. 

The  plots  for  U,V,W  and  P,Q,R  for  LOADEDREVERSAL  are  similar  to  those  for 
UNLOADEDREVERSAL,  except  for  the  value  of  W  during  the  reversal  (note:  they  are  plot¬ 
ted  on  different  timescales).  The  Euler-angle  plots  and  inertial  position  plots  are  notably 
different,  however.  As  is  shown  in  Figure  7.10,  the  LOADEDREVERSAL  involves  an 
increase  in  altitude  at  a  rate  of  about  30  feet/second  (as  compared  with  0  in  Figure  7.2)  and  a 
final  value  of  0  that  is  4  times  as  large.  Figure  7.11  shows  that  a  was  the  main  problem  for 
this  maneuver  (all  the  other  variables  tracked  their  commands  fairly  well).  The  plots  in  Figure 
7.12  show  that  Nz  did  remain  at  least  2  g  after  the  initial  bank  -  the  intersting  plot  is  the 
speed  change. 

The  cause  of  the  ramp  deceleration  can  be  seen  in  Figure  7.13,  where  it  is  shown  that  the 
throttle  turned  off  while  the  reversal  was  in  progress.  The  reason  for  the  throttle  shut  down  is 
revealed  in  Figure  7.16.  Note  the  plot  showing  the  actual  speed  squared  as  compared  with  the 
equilibrium  speed  squared  at  the  corresponding  a,p,9,<j>  points  along  the  trajectory.  The  prob 
lem  here  is  that  the  a,p,0,<j>  states  are  very  difficult  to  control  independently  near  a  =  0.  The 
throttle  saturated  trying  to  make  the  vehicle  slow  to  its  equilibrium  speed,  but  there  was  not 
nearly  enough  command  authority  to  reduce  the  speed. 

We  conclude  that  the  a.p.y.p  controller  was  not  a  good  choice  for  commanding 
LOADEDREVERSAL  —  the  failure  was  not  due  to  the  vehicle  or  its  flying  qualities.  It  is 
worth  pointing  out  that  the  a,p,0,0  controller  was  stable  at  every  point  along  the  maneuver,  as 
is  shown  by  the  plots  of  the  eigenvalues  (always  negative)  in  Figure  7.16.  This  is  a  case 
where  there  is  a  unique  physical  equilibrium  at  every  time,  but  control  saturation  prevented 
the  inversion  from  working  as  it  should. 

One  might  wonder  if  we  tried  to  bring  y  back  down  to  0  after  the  reversal  was  com¬ 
pleted.  The  answer  is  that  we  did  try  -  the  results  are  shown  in  Figures  7.17  through  7.24. 
Note  the  extreme  droop  in  N,  as  soon  as  the  command  for  y  decreased.  We  did  try  to  accom¬ 
modate  a  decrease  in  y  by  increasing  <J>,  as  can  be  seen  in  the  plot  for  p  and  MUCMD  in 
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Figure  7.19  but  it  did  not  work.  The  controller  was  simply  not  conditioned  well  enough  for 
this  land  of  maneuver. 
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Figure  7.13  Control  Activity  During  LOADEDREVERSAL 


Figure  7.14  A-Parameters  During  LOADEDREVERSAL 
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Figure  7.20  Accelerations,  Aero  Forces,  Alpha,  Beta  and  Speed  During  LOADED2REVERSAL 
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7.2  BARRELROLL 


The  maneuver  BARRELROLL  was  the  first  trajectory  generated  by  the  coordinated-flight 
U,P,Q,R  inversion  controller  described  in  section  6.  The  simulation  data  for  this  maneuver  are 
shown  in  Figures  7.25  through  7.32.  We  chose  the  command  profiles  for  U,  P,  and  Q  shown 
in  Figure  7.27.  The  strange  shapes  of  the  command  profiles  for  P  and  Q  after  t=12  seconds  is 
a  consequence  of  the  strategy  used  to  finish  the  maneuver  at  that  time  the  P  and  Q  commands 
were  used  to  drive  to  3<S0  degrees  and  6  back  to  its  equilibrium  value  at  t*3  seconds.  The 
commanded  R  was  computed  using  feedback,  approximately  by  the  formula  in  equation  6.4.2; 
the  term  associated  with  the  direct  side  force  due  to  rudder  was  neglected.  The  bandwidth 
chosen  for  R  was  only  one-tenth  the  bandwidth  chosen  for  P  in  this  maneuver.  In  fact,  this 
was  the  maneuver  which  led  to  the  bandwidth  parameters  discussed  in  section  6.4. 

The  simulated  responses  for  U,V,W,P,Q,R,<I>,8,'F,x,y,z  are  shown  in  Figure  7.25.  Note 
the  visible  V  component  that  arises  during  the  roll  in  the  U.V.W  plot  The  other  components 
of  the  speed  are  as  we  would  expect.  Outside  of  the  (small)  overshoots  in  the  angle-variables, 
the  other  plots  in  Figure  7.25  are  what  a  barrelroll  should  look  like. 

The  plots  in  Figure  7.26  show  the  change  in  altitude,  the  change  in  crossrange,  0,  and  ¥ 
as  a  function  of  time.  The  other  relevant  data  in  Figure  7.29  show  that  the  speed  remained 
nearly  constant  (as  it  was  supposed  to  do),  and  that  the  accelerations  were  all  tolerable.  The 
rather  large  (5  that  arose  did  cause  a  transient  Ny  of  about  0.1  g  in  magnitude.  That  could  have 
been  avoided  by  asking  for  a  higher  bandwidth  R  control  (relative  to  P). 

For  this  type  of  maneuver,  we  would  like  to  investigate  the  theoretical  implications  of 
equation  6.4.2  more  thoroughly  in  the  future. 
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Figures  7.33  through  7.38  depict  the  maneuver  DIVINGTURN,  which  we  chose  because 
of  its  dynamic  complexity.  A  ground  trace  of  the  maneuver  in  inertial  x-y  coordinates  is  shown 
in  Figure  7.33;  it  is  a  right-turn  (positive  y  is  to  the  right)  that  starts  at  roughly  x=2000  feet, 
y=0,  moves  out  along  the  positive  x  axis  and  then  turns  back  and  heads  into  the  negative  x 
halfspace.  The  change  in  aldtude  is  shown  in  Figure  7.34. 
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The  plots  of  the  vehicle  states  are  shown  in  Figure  7.35.  Note  that  U  drops  rapidly 
between  t=*4  seconds  and  t=7  seconds  —  a  change  of  more  than  300  feet  per  second  in  3 
seconds.  At  first,  we  might  suspect  that  the  pilot  experiences  some  very  unpleasant  accelera¬ 
tions  in  the  longitudinal  axis,  but  a  glance  at  the  N„  plot  in  Figure  7.36  shows  that  Nx  is 
almost  ruler  flat  at  -1.2  g  during  the  rum,  except  for  a  small  bump  between  6  and  8  seconds. 
Our  sign  convention  is  chosen  so  that  negative  N„  represents  a  force  pushing  the  pilot  back 
into  his  seat.  That  is  much  better  than  being  pulled  out  of  the  seat  at  3  g  ,  the  situation  easiest 
to  imagine  when  the  average  value  of  U  is  -3  g.  It  was  this  maneuver  that  alerted  us  to  the 
result  discussed  in  section  6.3,  that  the  aerodynamic  data  for  the  F-14  has  the  property  that  the 
aerodynamic  force  vector  is  almost  exactly  aligned  with  the  z-direction  of  the  aircraft  at  large 
a.  The  value  of  N,  shown  on  the  plot  is  almost  constant  at  -1.2  g  because  the  throttle  is 
saturated  at  full  power,  +1.2  g,  throughout  the  tum.  In  fact,  the  bump  between  t=6  seconds 
and  t=8  seconds  arose  because  a  exceeded  the  largest  recorded  datapoint  at  50  degrees  during 
that  time,  and  we  were  using  the  table  values  for  a  =  50  degrees. 

A  description  of  the  maneuver  is  as  follows:  the  pilot  begins  by  banking  the  aircraft  a  lit¬ 
tle  over  90  degrees,  then  (while  still  rolling)  he  pulls  back  hard  on  the  stick.  All  the  while  the 
throttle  is  set  at  full  power.  The  elevator  saturates  at  its  maximum  value  for  three  seconds,  but 
Q  quickly  turns  negative  because  the  inertial  terms  dominate  (see  the  software  user’s  manual, 
a  separate  document  prepared  under  this  contract).  As  Q  drops,  a  decreases  as  well  (note:  P 
is  0  during  the  time  when  Q  drops  off,  see  equation  6.4.4,  and  the  discussion  after)  until  P 
once  again  is  commanded  to  be  nonzero,  generating  (with  R)  some  extra  Q  to  make  Q 
increase  again  starting  at  t=8  seconds.  The  pitch  angle,  however,  decreases  steadily  to 
approximately  -60  degrees  until  time  t=  1 2  seconds.  Then  the  new  Q  command  and  the 
reverse  bank  command  at  time  t=12  seconds  brings  the  nose  of  aircraft  back  up  and  the  wings 
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Figure  7.38  Control  Activity  During  DIVINGTURN 


back  to  level.  Finally,  by  about  t=20  seconds,  the  aircraft  is  level  again  and  moving  at  the 
desired  speed. 


Conclusions 


From  the  three  types  of  maneuvers  flown  in  simulation  using  the  dynamic  inversion 
method,  we  were  led  to  several  of  the  candidate  nonlinear  flying  quality  parameters  discussed 
in  section  6,  and  we  were  forced  to  a  new  and  improved  dynamic  inversion  concept.  One 
parameter  provided  an  analytic  relationship  between  the  bandwidths  of  the  dynamic  P  and  R 
responses  with  a  restrictions  on  coordinated  flight,  another  allowed  us  to  understand  the  effect 
of  the  aerodynamic  force  vector  on  the  accelerations  experienced  by  the  pilot.  The  first 
maneuver  told  us  little  about  the  aircraft,  but  it  did  point  out  the  need  for  the  development  of 
the  controller  that  we  used  to  discover  the  two  important  relations  just  mentioned.  Other  ideas 
discussed  in  section  6  for  nonlinear  flying  qualities  also  came  out  of  analysis  of  these  simu¬ 
lated  maneuvers.  We  believe  the  dynamic  inversion  method  and  the  simulated  maneuvers  are 
a  very  useful  tool  for  developing  flying  quality  parameters. 


orxinuovwc*. 


W.sW.M 


IWKIV 


mw  htw mvmmncmnnncn  iivvuk 


References 


[Ang]  E  L.  Anglin,  "Static  Force  Tests  of  a  Model  of  a  Twin- Jet  Fighter  Airplane  for  Angles 
of  Attack  from  -10°  to  110°  and  Sideslip  Angles  from  -40 P  to  40°,"  NASA  TN  D-6425,  Aug. 
1971. 

[BCLA]  B.  Buchberger,  G.  E.  Collins,  R.  Loos,  R.  Albrecht,  "Computer  Algebra  Symbolic 
and  Algebraic  Computation,"  Springer- Verlag,  1983. 


[Cop]  W.  A.  Coppel,  "A  Survey  of  Quadratic  Systems,"  Journal  of  Differential  Equations  2, 
293-304,  1966. 

[KC]  I.  S.  Kukles  and  M.  Casanova,  "On  the  Distribution  of  Critical  Points  of  the  First  and 
Second  Group  (R). ,"  Izv.  Vyss.  Ucebn.  Zaved.  Mathematica,  No.  6,  88-97,  1964. 

[CM]  J.  V.  Carrol,  R.  K.  Mehra,  "Bifurcation  Analysis  of  Nonlinear  Aircraft  Dynamics,"  Jour¬ 
nal  of  Guidance,  Control,  and  Dynamics,  vol.  5,  no.  5,  Sept.-Oct.  ,  1982. 

[DMBS]  J.  Dongarra,  C.  Moler,  J.  Bunch,  and  G.  Stewart,  "LINPACK  Users’  Guide,  "  SIAM, 
Philadelphia,  1979. 

[El]  M.  R.  Elgersma  "Nonlinear  Control,  with  Application  to  High  Angle  of  Attack  and 
VSTOL  Right,”  Masters  Thesis,  University  of  Minnesota,  1986. 

[E2]  M.  R.  Elgersma,  "Control  of  Nonlinear  Systems,"  Ph.  D.  Thesis  ,  University  of  Min¬ 
nesota,  in  preparation. 

[E3]  B.  Etkin,  "Dynamics  of  Atmospheric  Flight,"  New  York:  Wiley,  1972. 

[G]  R.  Gilmore,  "Catastrophe  Theory  for  Scientists  and  Engineers,"  John  Wiley  and  Sons, 
1981. 

[HSM]  L.  R.  Hunt,  R.  Su,  and  G.  Meyer,  "Design  for  Multi-Input  Nonlinear  Systems," 
Differential  Geometric  Control  Theory,  edited  by  R.  Brockett,  R.  S.  Millman,  and  H.  J. 


Sussmann,  Birkhauser  1983  . 


[M]  J.  E.  Marsden,  "Elementary  Qassical  Analysis,”  San  Francisco:Freeman,  1974. 

[MKC]  R.  K.  Mehra,  W.  C.  Kessel,  J.  V.  Carrol,  "Global  Stability  and  Control  Analysis  of 
Aircraft  at  High  Angles  of  Attack,"  Cambridge:Scientific  Systems,  1977. 

[MMTJ]  D.  Mitchell,  T.  Meyers,  G.  Teper  and  D.  Johnston,  "Investigations  of  High-Angle- 
of-Attack  Maneuver-Limiting  Factors,  Part  3:  Appendices  -  Aerodynamic  Models,"  Technical 
Report  AFW AL-TR-80-3 141,  Part  3,  Systems  Technology,  Inc.,  December  1980. 

[PY]  M.  Pohst,  D.  Y.  Y.  Yun,  "On  Solving  Systems  of  Algebraic  Equations  via  Ideal  Bases 
and  Elimination  Theory,"  Proceedings  of  the  1981  ACM  Symposium  on  Symbolic  and  Alge¬ 
braic  Computation. 

[R]  A.  J.  Ross,  "A  Comparison  of  Analytical  Techniques  for  Predicting  Stability  Boundaries 
for  Some  Types  of  Aerodynamic  or  Cross-Coupling  Non-Linearities,"  article  17,  AGAARD 
CP-333,  1982. 

[SBDGIKM]  B.  T.  Smith,  J.  M.  Boyle,  J.  J.  Dongarra,  B.  S.  Garbow,  Y.  Ikebe,  V.  C.  Klema, 
C.  B.  Moler,  "Matrix  Eigensystem  Routines  -  EISPACK  Guide,"  New  York:Springer  Verlag, 
1976. 

[vdW]  B.  L.  van  der  Waerden,  "Algebra,"  Frederick  Ungar  Publishing  Co.  ,  1970. 

[W]  R.  J.  Walker,  "Algebraic  Curves,"  Springer- Verlag  ,  1978. 


[YSJ]  J.  W.  Young,  A.  A.  Schy,  and  K.  G.  Johnson,  "Pseudosteady-State  Analysis  of  Non¬ 
linear  Aircraft  Maneuvers,"  NASA  Technical  Paper  1758,  December  1980. 


APPENDIX  A 


TRAJECTORIES  IN  INERTIAL  SPACE 


The  equilibrium  velocities  will  give  trajectories  in  inertial  space  which  are  vertical  helices. 
One  way  to  show  this  is  as  follows. 

The  state  is  x  =  (V,p,a,<i>,0,'F,<I>,0)  .  At  equilibrium,  the  state  is  constant  so  two  of  the 
Euler  angles,  $  and  0  ,  are  constant  while  the  third  Euler  angle,  Y  ,  is  given  by 

'P(t)  =  'P(O)  +  4' t . 


The  transformation  matrix  between  body  axis  coordinates  and  earth  fixed  coordinates  is  given 
by  the  orthogonal  matrix  Upe<1)  =  Lq,  L©  where: 


Ly  - 


cosCP(t))  -sinOP(t))  O" 
sinOF(t))  cosOF(t))  0 
0  0  1 


cos(0)  0  sin(0) 
0  1  0 
-sin(0)  0  cos(0) 


1  0  0 


L*  = 


0  cos(G>)  -sin(O) 
0  sin(d>)  cos(<t>) 


The  velocity  vector  is  given  by 


u 

U 

V 

in  body  axis  coordinates,  and 

V 

w 

w 

VIPa 


where  1^  is  the  following  unit  vector 


1 


fia  ~ 


cos((3)cos(a) 

sin((J) 

.cos(P)sin(a) . 


So  in  earth  fixed  coordinates,  the  equilibrium  value  of  the  velocity  vector  is  given  by 


[i  a  cc 
VlpaJ=V  a  si 


cosCF(t))  -  b  sinOP(t)) 
sinCF(t))  +  b  cosCF(t)) 
c 


where 


b  I  =  Le 


The  right-hand  side  of  equation  A.1  is  usually  written  as 

cos(7)cosOJ'(t)-4'vb) 
V  cos(7)sinOI'(t)-'I'vb) 
-sin(Y) 

where  the  flight  path  angle,  y ,  is  given  by 


y  =  -sin 


r*(c)  =  cos-1  [Va2  +  b2] 


*  * 
%(,  =  sin-1  - 

V^Tb2 


(A.1) 


(A.2) 
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Note:  'Fyb  =  aircraft  heading  -  velocity  vector  heading. 


Integrating  (A.2)  gives  the  trajectory  in  the  earth  reference  frame. 


X«(t)  Xe(0)  sinCF(t)-^vb) 

Ye(t)  -  Ye(0)  =-^P-  -cosmt)-'Fvb) 
Ze(t)  Z*(0)  T  -tan(y)  'F  t 


This  defines  a  vertical  helix  of  radius  -qj-cos(y) 


sin(4/(0)-vt/vb) 

-cos(vF(0)-'Fvb) 

0 


(A.3) 
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APPENDIX  B 


DYNAMICAL  PROPERTIES  OF  FLIGHT  MANUEVERS 

GEORGE  R.  SELL* 

0.  Introduction.  In  this  report  we  ahow  that  the  atudy  of  the  dynamical  properties 
of  flight  manuevera  ia  governed  by  the  theory  of  time-varying  differential  equations.  The 
stability  properties  of  the  solutions  of  these  equations  are  determined  by  the  corresponding 
properties  of  the  limiting  eguotionj,  see  Sell  (1967b)  and  Sacker  and  Sell  (1977).  For  most 
flight  manuevera  these  limiting  equations  are  either  autonomous,  or  periodic  in  time.  In 
these  cases,  there  is  a  rid)  literature  for  describing  the  stability  properties.  A  brief  outline 
of  the  basic  theory  of  limiting  equations  is  induded 

One  of  the  first  objectives  ia  to  study  the  dynamical  properties  of  flight  manuevers  in 
the  vicinity  of  the  equilibrium  manifold.  The  Stable  Manuever  Theorem,  which  we  present 
here,  addresses  this  issue.  This  theorem  states  that  if  one  begins  a  flight  manuever  near  a 
strongly  stable  equilibrium  point,  and  if  the  manuever  input  is  close  to  a  nominal  input, 
then  the  aircraft  remains  near  the  equilibrium  manifold. 

The  report  contains  six  sections:  1.  Dynamic  Inversion,  2.  The  Equilibrium  Manifold, 
3.  Flight  Manuevera,  4.  Limiting  Equations,  5.  Applications  to  Flight  Manuevers,  and  6. 
Open  Problems. 

1.  Dynamic  Inversion.  Many  of  the  dynamical  properties  of  flight  manuevers  can 
be  understood  by  studying  a  model  of  the  dynamics  in  terms  of  an  ordinary  differential 
equation  with  a  control  parameter.  Typically  one  has  a  control-theoretic  problem  of  the 
form 

(1)  u)'  =  W(w,  u) 

where  u  €  RH  and  to  €  H*+".  We  assume  that  u  is  restricted  to  lie  in  a  fixed  open 
bounded  set  0  in  Rn.  Our  objective  is  to  describe  a  control  strategy  whereby  (1)  takes  on 
a  desired  form,  say 

w'  =  D(  is), 

where  D(w)  is  a  desired  vector  field  on  R"+"\  In  order  to  accomplish  this  we  need  to  solve 

(2)  W(w,u)mD(  is) 

for  a.  If  the  Jacobian  matrix  D»W(u>,u)  is  nonainguiar  and  if  D(w)  assumes  values  in  the 
range  *  (W(u>,u) :  u  €  0},  then  one  can  find  a  continuous  solution  u  *  u(u>)  of 

(2). 
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Since  the  atate  variable  w  ha a  higher  dimension  than  the  control  parameter  u,  the  above 
strategy  can  only  be  successful  if  some  restrictions  are  imposed  on  the  desired  function 
D(w).  This  means  that  some  of  the  variables  of  (1)  are  controllable  while  others  are  not. 
In  order  to  better  understand  this,  let  us  look  at  a  special,  but  rather  important,  case. 
Let  us  write  to  »  (x,y)T  where  a  €  A"  and  y  €  Rm.1  Then  (1)  takes  on  the  form 


r*'«/(*,y,«) 

l  »***(*.».  «o 


for  some  functions  /  and  §.  The  strategy  is  to  seek  to  control  the  x-variable  according  to 
the  desired  equation 

x'*D(x,y). 

In  order  to  do  this,  we  need  to  determine  a  control  strategy  u  by  solving  /(x,  y,  u)  »  D(zt  y) 
for  u.  As  before  this  leads  to  a  continuous  solution  u  «  U(x,  y)  whenever  the  Jacobian 
matrix  Dmf  is  nonsingular  and  £(x,y )  €  /(x,y,fl).  The  y-variables  are  not  controllable 
directly  by  this  strategy.  Instead  y  =  y(t)  must  be  a  solution  of  the  equation 


By  combining  this  we  see  that  (3)  takes  on  the  form 

(4)  *'  *  D(x,  y),  y'  «  y(x,y,[/(x,  y)). 


This  process  of  solving  for  u  is  referred  to  as  dynamic  inversion . 

The  dynamical  properties  of  uncontrollable  variables  y  are  determined  by  (4).  It  is 
these  properties  which  will  decide  whether  a  given  control  strategy  is  desirable,  and  we 
expect  that  the  same  properties  will  oftentimes  determine  the  flying  qualities  of  an  aircraft 
and  a  collection  of  flight  manuevers.  We  will  illustrate  this  in  a  moment,  after  we  define 
the  Equilibrium  Manifold  for  (1). 

2.  The  Equilibrium  Manifold.  In  order  to  simplify  our  treatment  a  bit,  we  will 
assume  a  special  form  of  (1)  which  is  given  by 


(5)  w'  *  F(u>)  +  G(w)H(w,u), 

where  the  functions  F,  G  and  H  are  smooth  functions  with 

F  :Rm+m  -  Rn+m 
G  :Rm+m  -+  L(Rm,Rn+m) 

H  :Rm+m  x  Rn  —  Rn 

1  Mon  |«wnll)f  on*  could  introduce  local  coordinates  w  c  (x,  y)T  and  rostrict  *  to  li«  on  some  given 
manifold  ht  in  Tho  function  D  is  than  a  daoired  vector  geld  on  J it,  and  y  denote*  the  normal 
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aud  £,(A",  An+’")  denote*  the  space  of  linear  mappings  from  Rn  into  A"+m.  In  other 
words,  for  each  w  €  Rn+m  and  u  €  A",  F(u>)  is  a  (n  +  m)  x  l  column  vector,  H (u>,  u)  is 
a  n  x  1  column  vector,  and  G(w)  iia(n  +  m)xn  matrix.3 

Let  ft  denote  the  subset  of  A"*"*  where  the  rank  of  the  Jacobian  matrix  DwG(w)  is 
n,  and  assume  that  0  is  a  nonempty  open  set  A".  This  permits  one  to  introduce  local 
coordinates  w  m  (x,y)T  so  that  (5)  takes  on  the  equivalent  form 


f  *'  *  /i(*,»)  +  fi(*.  V)H(*>  V.«0 

\  v'  *  /*(*.v)  +  ya(*,y)ff(*,y.tO- 


where  g » is  an  n  x  n  matrix,  and  ya  is  an  2  x  n  matrix.  We  assume  that  the  local  coordinates 
have  been  chosen  so  that  gi  is  invertible.3  Here  x  denotes  the  controllable  variables  and  y 
the  uncontrollable  variables.  Next  we  define  a  manifold  Mo  for  (6)  as 


Mo  *  {(*,V)  €  A"*"' :  /i(x,y)  +  yi(*,y)A(*,y,u)  «=  0}. 


The  manifold  Mo  is  invariant  for  (6),  and  on  Mq  one  has 

(7)  H  (x,  y,  ti)  =  -yf  *(*.  v)/i(*.  y)- 
By  inserting  (7)  into  (6)  we  obtain 

(8)  y'-AXx.y) 


where 

A'(x,y)  =  f7{z,v)  -  97(x,y)gi'{x,v)fi(x,y) 
and  x  =  constant.  On  the  invariant  manifold  Mo  equation  (6)  reduces  to 


(9) 


(  v'  •  A'(*,y)  *  Mx,y)  -  92(x,y)g^l(xty)fl(x,y). 


By  fixing  H  so  that  (7)  holds,  the  controllable  variable  x  is  held  constant.  The  function 
H  may  change  with  time,  but  this  only  reflects  the  adjustment  in  H  caused  by  the  time 
variation  of  the  uncontrollable  variable  y,  which  is  a  solution  of  (8). 

The  Equilibrium  Set  for  (6)  is  definded  as  the  set 


E  *  {(*.y)  €  Mo  :  /a(*,y)  +9t(*,v)R  «  0), 

*We  «iU  esnesntrsts  our  attention  on  the  special  model  equation  (5).  The  general  dynamical  theory 
wa  describe  here  is  not  dependent  on  the  fact  that  the  control  in  (5)  enters  as  a  linear  factor  H  .The  only 
role  that  the  special  form  (S)  plays  in  this  report  is  to  simplify  some  of  the  algebraic  considerations  which 
arise  below. 

*This  is  possible  because  of  the  rank  condition  on  the  Jacobian  matrix  DVC. 


*• 

Figure  1.  The  Equilibrium  Manifold  near  a  point  x0. 
where  If  ia  given  by  (7).  Notice  that  (x,y)  €  E  if  and  only  if  (x,y)  €  Mo  and 

(10)  K(*,V)  «  /a(*.v)  “  9i(*'V)9il(*>v)fi(*,v)  =  0. 

If  the  Jacobian  matrix  DtK  ia  nonaingular,  then  one  can  aolve  (10)  locally  for  y  =  e(x) 
where  c(x)  is  continuoua  in  x.4  It  is  possible  that  for  a  given  x,  equation  (10)  may  have 
several  solutions  y  «  e(x),  see  Figure  1. 

The  Equilibrium  Manifold  is  a  subset  of  the  equilibrium  set  and  is  defined  by 

M  =  {(x,y)  €  E  :  DtK{x%y)  is  nonsingular}. 

For  each  (io, Vo)  €  M  there  ia  a  neighborhood  U  of  x0  and  a  C1 -function  e  :  U  — *  Rm 
such  that  (x,e(x))  €  M  for  x  €  U  and  e(xo)  *  yo-  This  means  that  in  the  vicinity  of 
every  point  on  the  equilibrium  manifold  M,  the  equilibrium  set  E  agrees  with  M,  and  E  is 
locally  a  smooth  manifold.  The  behavior  of  E  near  points  (x,y)  €  E  where  DtK(x,y)  is 
singular,  can  be  very  complicated.  These  singular  points,  which  can  be  bifurcation  points, 
will  not  be  analyzed  in  this  report. 

Since  the  equilbrium  points  for  (6)  are  generally  not  isolated,  it  is  not  possible  for  them 
to  be  asymptotically  stable  in  the  full  dynamics  (6).  However  it  is  possible  that  equilibrium 
points  of  the  y-  equation  (8)  are  asymptotically  stable.  Because  of  the  importance  of  this 
fact,  we  introduce  the  following  stability  concept:  We  shall  say  that  at  point  (xo,y0)  is 
strongly  stable  if  (xo,yo)  is  an  equilibrium  point  for  (6),  and  yo  is  asymptotically  stable 
for  equation  (8). 

Assume  for  the  moment  that  for  each  x  €  Rn,  equation  (8)  has  an  isolated  equilibrium 
point  c(x)  and  that  e(x)  varies  continuously  in  x.  Assume  further  that  over  some  interval  of 
time  I  the  controller  u(t)  is  chosen  so  as  to  maintain  a  given  value  for  x  and  that  y  =  e(x) 
for  £  €  /•  (Since  x  and  y  are  constant  on  /,  it  follows  from  (7)  that  H  is  constant  as  well.) 

4Aa  argument  for  showing  the  exist  sites  of  e  solution  c(s)  is  given  below  as  a  part  of  the  proof  of  the 
Stable  Manuever  Theorem. 


Next  we  ask,  what  would  happen  if  a  email  disturbance  is  introduced  into  this  system? 
Ciparly  this  will  be  strongly  influenced  by  the  stability  properties  of  the  equilibrium  point 
e(x).  If  the  linearization  of  the  vector  field  at  e(x)  gives  rise  to  an  eigenvalue  with  positive 
real  part,  then  one  can  expect  an  abrupt  change  in  the  response.  One  would  like  (0)  to 
have  strong  stability  properties  and  to  be  free  of  destabilizing  bifurcations.  We  shall  look 
into  this  further  in  Section  5. 

Fbr  equation  (6)  it  is  convenient  to  separate  the  problem  of  controllability  into  two 
parts.  Since  a  enters  this  problem  through  the  function  H  only,  one  can  think  of  H  itself 
as  a  controller.  Fbr  example,  the  aircraft  pilot  may  select  a  specific  value  H0  for  H  at  a 
given  time  t,  then  the  onboard  computer  would  solve  the  equation 

(11)  lf(xty,u)-lfo 


for  u  to  drive  the  rudder,  thrust,  etc.  The  solvability  of  the  last  equation  for  u  depends,  of 
course,  on  the  Implicit  Function  Theorem,  and  usually  requires  the  Jacobian  matrix  DUH 
to  be  nonsingular.  We  emphasize  that  the  dynamical  properties  we  are  studying  here  are 
independent  of  whether  or  not  (11)  is  solvable  for  u. 

So  for  we  have  been  tacitly  assuming  that  any  point  in  the  (x,  y)- spare  is  attainable 
via  a  flight  manuever.  This  is  by  no  means  the  case.  Furthermore,  the  coordinates  for  the 
problem  need  not  be  Euclidean  coordinates.  Angular  variables  are  natural  in  many  cases. 
What  this  means  is  that  one  should  consider  the  original  problem,  either  (1)  ro  (5),  to  be 
given  on  a  prescribed  subset  A  of  some  manifold  M.  We  shall  refer  to  A  as  the  attainable 
set.8 

In  applications  developed  by  others  on  this  project,  one  has  n  =  4,m  =  2,  and  conse¬ 
quently  (8)  represents  an  ordinary  differential  equation  in  the  plane  R 7 .  The  equilibrium 
manifold  M  in  this  case  has  dimension  4. 

3.  Flight  Manuevers.  We  define  a  flight  manuever  to  be  the  response  of  an 
aircraft  to  a  time-varying  controller  on  some  time  interval  /.  For  the  model  equation 


,  12)  f  *  /»(*.  V)  +  f  »(*,  V)H 

\v'*/a  (*,v)  +  Pa(x,y)H 

where  H  “  Jf(t)  *  if(x,y,ti,<)  is  now  a  time-varying  input,  a  flight  manuever  is  a  solution 
(*(0»lK0)  (12),  for  t  €  /  ■*  (a,  i).  The  function  H  is  referred  to  as  the  manuever 

input.  We  will  assume  that  the  inputs  H(t)  are  piece-wise  continuous  functions  of  t  on  /. 
Furthermore,  we  restrict  our  attention  to  inputs  H{t)  for  which  the  response  (x(f),  y(f)) 
remains  in  the  attainable  set  A. 

*  Event  hough  vt  shall  eootiaus  to  writ*  our  aquations  in  terms  of  Euclidean  coordinates,  this  is  really 
not  essential  tor  our  theory. 
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There  are  two  manuevers  which  are  of  special  interest  for  flight  control.  The  first  of 
these  is  s  basic  manuever  on  an  interval  /.  This  manuever  occurs  when  there  exists 
* a,ya  such  that 

*(0  -  -»,•'(*..  »(«))/.(*..»«))  <  €  /, 

where  y(t)  is  the  solution  of  (8)  with  y(e)  *  ya.  In  a  basic  manuever  the  quantity  xa  is 
referred  to  as  the  base  state  of  the  aircraft.  For  a  basic  manuever  the  system  (12)  reduces 
to  the  system  (9).  These  manuevers  are  quite  common  for  aircraft  motion.  They  occur, 
for  example,  when  the  aircraft  is  ascending  (or  descending)  with  a  fixed  angle  of  attack, 
or  when  the  aircraft  is  turning  with  the  controls  held  fixed. 

The  second  manuever,  which  we  call  an  advanced  manuever,  occurs  when  H(t)  is 
continuous  at(»s,i  and  one  has 

m  -  -9r,(*(<),v(0)/i(*(0.v(0)  * - 

where  (x(t),  y(t))  is  a  solution  of  (12)  on  /.  Advanced  manuevers  occur,  for  example,  during 
a  barrel  roll,  or  when  the  pilot  is  changing  the  thrust  vector.  Any  aircraft  trip,  from  takeoff 
to  landing,  is  a  sequence  of  basic  manuevers  interspersed  with  advanced  manuevers. 

There  is  a  convenient  way  to  distinguish  between  basic  and  (non- basic)  manuevers. 
Asstime  that  the  manuever  input  H  satisfies 

(i3)  m  -  »«)/.(*<<).  »(<».  <  e  i- 

It  then  follows  that  the  x -equation  reduces  to  x1  —  0,  or  x(<)  =  xa,  a  constant,  for  t  €  I. 
In  other  words  a  manuever  is  a  basic  manuever  if  and  only  if  (13)  is  satisfied.  The  extent 
to  which  H(t)  deviates  from  satisfying  (13)  is  a  measure  of  the  strength  of  a  manuever. 
We  can  quantify  this  by  defining  the  norm,  or  strength,  of  a  manuever  N(H)  by 

W.H)  -  .up  f|ff(!)  +  *'<*(<).  v(0)/.(*(0.»(<))l  :<€/)• 

Notice  that  one  has  N(H)  =  0  if  and  only  if  H  is  the  input  for  a  basic  manuever.  In  the 
Stable  Manuever  Theorem,  which  we  give  below,  we  show  that  if  a  manuever  begins  near 
a  strongly  stable  equilibrium  point  on  the  equilibrium  manifold  M  and  if  N(H)  is  small, 
then  the  aircraft  remains  near  the  strongly  stable  equilibrium  points  on  M  throughout  the 
entire  manuever. 

An  advanced  manuever  can  be  viewed  as  trajectory  which  changes  the  base  state  of  an 
aircraft  from  xa  to  x».  The  point  to  emphasize  is  that  the  full  equations  (12)  are  needed 
to  describe  the  dynamical  behavior  of  an  aircraft  during  an  advanced  manuever.  However 
at  time  t  *  a,  b  these  equations  reduce  to  (9).  Since  an  advanced  manuever  ends  with  (9) 
being  satisfied  at  t  m  &,  it  follows  that  the  equations  (12)  are  asymptotically  autonomous. 
We  shall  present  next  a  brief  review  of  limiting  equations  and  asymptotically  autonomous 
equations  in  the  next  section. 
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4.  Limiting  Equations.  We  begin  with  a  nonlinear  time-varying  differential  equa¬ 
tion  x'  «  /(x,f),  where  /  €  C(RN  x  R,RN)  sad  C(RN  x  R,RN )  denotes  the  space  of 
continuous  functions  from  RN  x  R  to  RN.  In  order  to  simplify  the  discussion  we  will 
assume  that  all  functions  /(*,<)  considered  in  tills  section  arc  Lipscliitz  continuous  in  x, 
and  that  every  solution  of  x'  *  /(x,f)  is  defined  for  all  t  € 

For  /  €  C(RN  x  RtRN)  we  define  the  translation  fr  by  /r(x,t)  *  /(x,t  +  T)-  The 
mapping 

e(f,r)  =  fr 

defines  a  flow  on  C(Jlw  x  R,RN).  If  we  let  d(*i/<0  denote  the  solution  of  x*  —  /(x,t) 
satisfying  d(*./»0)  “  x,  then 


*(*./.*)«(*(*./.  r)./r) 

is  a  (skew-product)  flow  on  RN  x  C(/?v  x  R,  RN). 

In  the  flow  a  on  C(RN  x  R,  RN )  we  define  the  hull7  of  /  by 

*(/)-Cl  {/r  :r€R) 


and  the  positive  hull  by 

H+U)  «  Cl  (fr  :  r  >  0}. 

The  w-limit  set  of  a  function  /  €  C(RN  x  R ,  #*)  is  given  by 

«(/)  =  nr>„^+(/r). 

The  limiting  equations  for  /  €  C(RN  xR,RN)  is  defined  as  the  collection  of  all  equations 
x'  wt  g(x,  t)  with  g  €  ft(/).  The  following  theorem  gives  a  useful  sufficient  condition  for  the 
collection  of  limiting  equations  to  be  nonempty  and  compact.  The  proof  of  this  theorem 
is  given  in  Sell  (1967ab)  and  is  based  on  the  Ascoli-Arzela  Theorem. 

Limiting  Equation  Theorem.  Let  f  e  C(RN  x  R,RN)  be  such  that  for  every 
compact  set  K  C  RN  there  is  a  constant  k>0  and  a  function  6  =  6(e)  with  the  following 
two  properties: 

(1)  /(x,t)  is  Lipsehits  continuous  in  x,  uniformly  in  t,  i.e.  |/(x,t)  -  /(y,t)|  <  h|x  -  y| 
tor  nil  x,y  €  /C,t  €  it 

(2)  /(x,t)  is  uniformly  continuous  on/f  xA,  i.e.  |/(x,s)  —  /(y,t)|  <  c  for  nil  x,y  € 
A',s,<  €  fl  with  |x  -  y|  <  6(e)  and  |t  -  s|  <  6(e). 

*The  assumptions  on  the  LipachiU  continuity  of  /  and  the  global  existence  of  solutions  of  r'  as  /(*,  <) 
can  be  dropped.  See  Sell  (1997ab,1973),  MUier  and  Sell  (1970)  and  Seeker  and  Sell  (1977)  for  details. 

7 The  closure  operation  Cl  used  here  is  the  closure  in  the  topology  of  uniform  convegcnce  on  compact 
subsets  of  RH  x  R. 
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Then  Cl(f)  is  a  nonempty,  coin  pert,  connected  subset  of  C(RN  x  R,RN).  Furthermore 
every  g  €  fl(/)  satisfies  the  Lipscbitx  condition  (1)  above. 

The  maiu  advantage  of  the  limiting  equations  in  Unit  one  out  use  tltn  standard  theory 
of  u^limit  Mti  to  describe  the  behavior  of  solutions  of  the  origiual  differential  equation 
x*  ■»  /(*,<)  os  t  — *  oo.  In  particular,  let  /  €  C(RN  x  R,RN)  and  let  ^(*,/,t)  be  a 
solution  of  x'  «  /(x,t)  that  stays  in  a  given  compact  set  /f  C  /?N  for  t  >  0.  Assume 
further  that  /  satisfies  the  hypotheses  of  the  Limiting  Equation  Theorem.  Then  the 
motion  x(x,/,r)  remains  in  compact  set  in  RN  x  C(RN  x  R,RN).  Furthermore  for  every 
sequence  r„  — »  oo  there  is  a  subsequence,  which  we  denote  again  by  r„,  and  a  point 
(y,g)  €  RN  x  C(RN  x  A,/**)  such  that 

*(*./.  r„)  —  V 
/r.  —  9 

4{*,f,t  +  rn)  -*  t(y,g,t) 

where  the  last  limit  is  uniform  on  compact  subsets  of  R. 

There  are  two  special  situations  concerning  limiting  equations,  which  arise  in  the  theory 
of  flight  manuevers.  We  say  that  a  function  /  €  C(RN  x  R,RN)  is  asymptotically 
autonomous  if  /  satisfies  the  hypotheses  of  the  Limiting  Equation  Theorem  and  if  the  u- 
limit  set  of  /  consists  of  one  point,  say  fl(/)  ■=  (g).  Since  the  u-limit  set  is  invariant  under 
the  flow  a  on  C(RN  x  R,RN),  it  follows  that  gT  =  g  for  all  r  €  /?,  i.e.  y  is  autonomous 
(independent  of  time).  Similarly  a  function  /  6  C(RN  x  R,RN)  is  asymptotically 
periodic  if  /  satisfies  the  hypotheses  of  the  Limiting  Equation  Theorem  and  if  ft(/) 
consists  of  a  single  periodic  orbit.  In  this  case  g  €  f !(/)  is  periodic  in  t. 

Let  x'  /(x,f)  be  an  asymptotically  autonomous  differential  equation  with  limiting 
equation  x*  *  y(x).  Let  xt  be  an  asymptotically  stable  equilibrium  point  for  x'  =  y(x). 
Then  one  can  show  that  there  is  a  neighborhood  U  of  xo  such  that  for  all  *i  €  U  the 
solution  ^(*i ,/,<)  of  x'  *  /(x,f)  satisfies  ^(xj,/,f)  — *  x0  as  t  — »  oo,  see  Markus  (1956). 
In  other  words,  the  solutions  of  x*  «  /(x  ,t)  beginning  in  U  are  stable,  as  a  matter  of 
fact,  they  are  asymptotically  stable.  The  stablility  properties  of  asymptotically  periodic 
equations  are  similar,  see  LaSalle  (1962)  and  Sell  (1966). 

In  this  report  we  shall  restrict  our  study  to  flight  manuevers  which  are  asymptotically 
autonomous.  Eventhough  it  may  appear  that  all  flight  manuevers  are  asymptotically 
autonomous,  this  need  not  be  the  case.  For  example,  if  one  wishes  to  study  the  effects  of 
small  random  disturbances  (i.e.  noise)  on  the  aircraft,  then  one  may  need  the  full  theory 
of  limiting  equations.  The  references  cited  below  offer  a  good  introduction  to  this  theory. 

For  applications  to  the  study  of  the  dynamical  properties  of  flight  manuevers,  one 
should  note  that  the  limiting  behavior  of  a  function  /  may  be  assumed  in  finite  time.  In 
particular  if  /  €  C(RN  x  R,RN)  has  the  property  that  there  is  an  autonomous  function 
g  €  C(RN  x  R,RN)  and  a  time  T  such  that  f(x,t)  *  g(x)  for  all  t  >  T,  then  /  is 


asymptotically  autonomous  and  fl(/)  =  {y}.  Similarly  if  /  €  C(RS  x  R,RN)  has  the 
property  that  there  is  a  time-periodic  function  g  €  C(RN  x  Ry  RN)  and  a  time  T  such  that 
/(*.0  *«(*.<)  fc*  all  <>T,  then  /  is  tutymptoticnlly  periodic  and 


ft(/)=(Sr:0<T<P} 


where  P  is  a  time-period  of  g ,  see  Sell  (10C7ab). 

5.  Applications  to  Flight  Manuevers.  During  an  advanced  manuever  on  an  in¬ 
terval  /  *  (a,  b)  the  equations  of  motion  are  given  by 


f  *'  =  /i(*.v)  +  $i(*.v)ff 
1  v'  =  /j(*.y)  +  sa(*,v)ff- 


The  motion  begins  at  (xa,ya)  at  time  t  «s  a  and  ends  at  (x4,y4)  at  time  <  *=  i.  Equation 
(12)  is  an  asymptotically  autonomous  system  with  the  limiting  equation  given  by 


(9) 


x'  =  0 

y'  as  K(xt,  y)  =  /2(x4,  y)  -  ya(x4,  y)gf '(xs,  y)/i(x4,  y). 


As  a  matter  of  fact  (12)  reduces  to  (9)  for  t  >  b.  Clearly  the  location  of  y4,  the  value  of  y 
at  time  t  =  i,  can  play  a  significant  role  on  the  dynamics  of  the  manuever.  Let  us  look  at 
a  few  examples. 

Example  1:  Assume  that  y4  lies  in  the  basin  of  attraction  of  an  asymptotically  stable 
equilibrium  point  e(x4)  of  y'  =  A"(x4,  y).  In  this  case,  y(t)  is  attracted  towards  e(x4)  for 
t  >  b.  If  ya  is  close  enough  to  e(x*),  then  the  transient  behavior  will  be  short-lived,  and  for 
all  practical  purposes,  one  would  have  y(<)  <-  e(x4)  after  some  short  time  interval.  (This 
situation  is  quite  common  and  very  likely  is  the  outcome  of  most  flight  manuevers.)  On 
the  other  hand,  if  y4  is  far  away  from  e(x4),  then  the  transient  behavior  will  persist  for  a 
long  time,  and  the  aircraft  could  be  under  the  influence  of  the  transient  part  of  y(f)  when 
the  next  advanced  manuever  commences.  For  example,  if  y*  is  close  to  the  boundary  of  the 
basin  of  attraction,  then  one  would  expect  that  the  dynamical  behavior  of  the  boundary 
set  will  have  greater  influence  on  the  aircraft  than  that  of  the  equilibrium  point  e(x4). 

Example  2:  Assume  that  y4  does  not  lie  in  the  basin  of  attraction  of  any  asymptotically 
stable  equilibrium  point.  For  instance,  y4  may  lie  on  a  periodic  orbit,  or  in  an  unstable  or 
chaotic  portion  of  the  y- space.  This  could  have  very  serious  consequences  for  the  aircraft. 
One  would  not  want  to  implement  such  a  manuever  without  detailed  knowledge  of  the 
underlying  dynamical  properties  of  (9). 

Because  of  the  considerations  raised  in  the  last  two  paragraphs,  one  of  the  goals  in 
the  study  of  the  dynamics  of  flight  maneuevers  is  to  understand  the  dynamical  properties 
of  y'  «  K(n,y)  in  every  fiber  A(x4),  where  A(x)  is  defined  to  be  those  y  for  which 


(x,y)  €  A.  The  set  A(z)  is  simply  the  fiber  of  the  attainable  set  A  over  the  point  x. 
Because  of  our  assumption  that  every  input  H(t)  yield  an  attainable  response,  one  must 
have  y(<)  €  /t(x(t)),  for  all  t  €  /.  One  would  like  to  show  that  A(x*)  has  additional 
dynamical  properties.  For  example,  one  may  want  to  show  that  if  an  advanced  manuever 
ends  in  a  position  (x*,  y»),  then  the  trajectory  y(t)  of  y'  K(xt,  y)  with  y(b)  =  y»  remains 
in  A(x»)  for  all  t>  b. 

While  the  overall  dynamical  properties  of  flight  manuever*  can  be  very  complicated, 
there  is  one  practical  situation  which  is  amenable  to  analysis.  In  this  case  we  consider 
an  advanced  manuever  which  begins  at  a  point  (xai  ya)  which  is  close  to  the  equilibrium 
manifold  A#.  This  means  that  ya  is  close  to  some  equilibrium  state  e(xa)  of  y*  *■  J?(xa,  y). 
We  assume  that  e(xa)  is  asymptotically  stable  and  that  the  manuever  input  H(t)  is  close 
to  the  nominal  value  — yj’l(x(t),y(l))/i(x(f).y(<))*  i.e.  the  norm  Af(Jf)  i*  small.  We  will 
show  that  the  terminal  value  (x*,  y»)  is  close  to  the  equilibrium  manifold  Af  and  that  y»  is 
close  to  an  asymptotically  stable  equilibrium  point  e(x»)  of  y'  =  A'(x*,  y).  More  precisely 
we  will  prove  the  following 

STABLE  Manuever  Theorem.  Let  (xa,e(xa))  be  a  strongly  stable  equilibrium  point 
on  the  equilibrium  manifold  M,  i.e.  e(xa)  is  an  asymptotically  stable  equilibrium  point  for 
y‘  wt  K(xm, y).  Then  there  exists  positive  constants  <c,ci,ca,ci  and  e0  such  that  for  any 
«,0  <  e  <  «o»  the  following  holds:  Let  H(t)  be  the  input  for  any  manuever  on  any  interval 
I  m  (a,  b)  where  the  following  conditions  are  satisfied; 

(1)  The  manuever  begins  at  (xc,ys)  where  ]ya  —  e(xa)|  <  cj«. 

(2)  One  has  |x(t)  —  xa|  <  ej<  for  a  <  t  <  b. 

(3)  One  has  |Jf(l)  +  pr>(x(0iV(())/i(s(0*y(t))l  ^  N{H)<cst  fora<t<b. 

Then  for  each  r  €  /  there  exists  a  strongly  stable  equilibrium  point  (x(r),  e(x(r)))  of  (12) 
such  that  |y(r)  —  e(x(r))|  <  Kt  for  all  r  €  I. 

Before  proving  this  result  one  should  note  that  we  do  not  assume  anything  concerning 
the  length  of  the  time  interval  I  =  (a,  b).  The  interval  length  can  be  large.  The  constants 
*,ci,ca,cs  and  «o  are  independent  of  the  interval  I.  As  a  result,  this  theorem  does  apply 
to  any  sequences  of  basic  and  advanced  manuevers  which  satisfy  (2)  and  (3). 

Proof  of  Stable  Manuever  Theorem.  There  is  no  loss  in  generality  in  assuming  that 
(x*,e(xa))  w  (0,0).  (Indeed  if  this  were  not  the  ease,  then  the  change  of  variables 

*-»x  +  xa 

V-»V +  «(*•) 

would  result  in  a  new  system  of  differential  equations  with  the  desired  property.)  As  a 
result  one  has  A'(0,0)  *  0  and  the  matrix  A  =  DtK( 0,0)  is  stable.  The  latter  means  that 
there  exist  K\  >  0  and  A  >  0  such  that 

(14)  |e*‘|  <  /f,e"Al,  t  >  0 
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Next  define 

A(0  -  £  (t)  +  g;\x(t), y(t))/i(x(<),  y(f )),  t  €  /• 

Then  (12)  reduce*  to 

fx'(t)«y,(x(f).y(t))A(t), 

1  VV)  «  A'(*(0.lK0)  +  *(*(0.  v(0)A(t). 

Furthermore,  for  <  €  /  one  ha*  |A(<)|  <  N(H),  the  norm  of  the  manuever. 

Our  first  step  is  to  show  that  there  exist*  a  continuous  solution  y  =  e(x(t))  of  the 
equation 

A'(x(<),y)-0,  tel. 

This  is,  in  fact,  a  consequence  of  the  Implicit  Function  Theorem.  Since  we  also  need  an 
estimate  of  the  size  of  e(i(t))  we  will  present  the  details  here.  Define  £(x,y)  by 

£(*•  v)  “  y  —  X”1  A'(x,y). 

Then  £y£(0,0)  *  0.  Next  fix  ao  >  0  and  bo  >  0  so  that 

(15)  i/)|  <  j.  |x|  <  ao,  Ivl  <  V 

Let  M\ ,  A/j ,  A/j ,  Mi  be  constants  so  that 

|£(x,0)  —  £(0,0)|  <  A/j|x|, 

a6)  \D,K-\z,y)DtK{z,y)\  <  Ma, 

M*.v)l  <  Mi, 

.  M*,y)l  <  Mit 

for  |x|  <  ao,  |y|  <  V  By  making  ao  smaller,  if  necessary,  we  can  assume  that  2Afiao  <  f»o  - 
Because  of  (15)  the  mapping  y  — »  £(x,y)  is  a  strict  contraction,  and  consequently  the 
equation  y  *  £(x,y)  ha*  a  fixed  point  y  ■=  e(x)  for  |x|  <  to-  Furthermore  y  **  e(x)  is 
also  a  solution  of  J?(x,y)  m  0,  and  e(x)  is  a  CJ -function  of  x  which  satisfies  the  following 
estimates:  With  y  *  e(x)  *  £(x, y)  and  0  *  £(0,0)  one  ha* 

|e(x)|  «  |£(x,y)  -  £(0,0)  +  £(x,0)  -  £(x,0)| 

<  |£(x,y)  -  £(x,0)|  +  |£(x,0)  -  £(0,0)| 

<  jl«(*)l  +  A/i|x|. 

Hence  |e(x)|  <  2M\\z\  <  2M\Oo  <  b0.  Consequently  for  x  *  x(t)  one  has 

|e(x(t))|<2M1|x(t)|<2M1ca« 
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provided  «<«©<§*•  Furthermore  one  has 

je(*(0)  -  D,'U(t)),\t) 

m  DMe(x(t))gi(x(t)ty(t))A(t). 

If  one  differentiate*  e(x)  «  e(z)  —  A_,A'(x,e(x))  with  respect  to  x,  one  obtain* 

D, tf(x,e(x))D.e(x)  *  -I>.tf(x,e(x)) 

and  consequently  from  (16) 

|£>«e(x)|  <  |DtA'-*(xte(x))D,tf(x,e(x))|  <  M7. 

Next  we  define  x(<)  »  y(t)  —  e(x(t)),  and  for  the  remainder  of  the  proof  we  let  x  » 
x(t),x  »  x(t),y  =  y(t)  and  e  *  e(x(t))-  Then 

(17)  x' *s  Ax  -  I(x,  x), 

where  L  is  given  by 

L(x,  x)  *  A(x,  y)  -  A*  +  *3(x,  y)A  -  e' 

sa  A(x,  y)  -  Ax  +  yj(x,  y)A  -  I>,e(x(t))y,(x,  y)A. 


Now  fix  17  so  that 

(18)  . ,^2, 

where  Ki  and  A  are  given  by  (14).  Next  choose  ai,6j,0  <  oj  <  a0,0  <  6j  <  6o  so  that 


(1») 


p,jr(*,i()-o,^(0,0)|<7 


for  |x|  <  Si.lvl  -  hi*  Since  K(x,e(x))  *  K( 0,0)  =*  0  we  obtain  the  following  from  the 
Taylor  expenaian: 

K(x%  y)  -  Ax  *  K(x,  y)  -  A'(x,  e(x))  -  Ax 

«  A'(0,0)  +  (j  D,K(x,e  +  9x)d9Sj  x  -  Ax 

,K(x,t  +  6x)  -  D,*(0,0)l<tf)  x. 
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Assume  that  «o  satisfies 


<o  <  min 


If  1*1  <  <o  we  have  |e|  <  2A/j«0,  and  therefore  if  |r|  <  ^-  we  have  |e  +  0z|  <  b\.  FVom  (19) 
we  get 


(20)  !*(*,») 

for  |z|  <  e0,|y|  <  fcj.  FVom  (16)  we  obtain 

l*i(*,v)A|  <  A/aJV(ff)  <  A/jCj< 

1  '  \Dge(x)gl(x,v)A\<M2M4N(H)<MiM4ete. 

Now  (20)  and  (21)  imply  that 


(22)  |I(x,*)|<^|z|  +  M5e 

for  |*|  <  «  <  eo,|x|  <  where  A/a  »  A/jCj  +  Mj M4cj.  By  the  Variation  of  Constants 
Formula,  the  solution  of  (17)  is  given  by 

r(0  =  e*(,-*>*.  +  J*  eA^L(x(s)tx(3))ds. 

This  means  that  x(t)  is  a  fixed  point  of  the  operator 

F(z)(t)  =  e>l<‘->*.  +  j\A^-$)L(x(s)tx(s))ds. 

Assume  that  |z(j)J  <  for  a  €  /•  Then  by  using  (14)  and  (22)  we  have 

|f(«)(<)l  <  *■.«'»<— >|».|  +  /  «-*"->(«,.  +  !||j(j)|)lt> 

<  -  .-»<■-))  +  k.t, J 

Since  |z«|  =  |y«  -  e(x«)|  <  Cj«  <  CiC0  we  get 


(23) 


|f(«X0l  S  u,t  +  k,t)  f  «-*<— >  |,(.)|* 


where  A/«  *  K\C\  +  and  ||z||  *  sup  {|z(z)| :  s  €  /).  Since  ))r||  <  it  follows 

from  (18)  that  |F(zX<)l  -  In  order  to  assure  that  |F(zX<)l  S  V  **  re<lu^re  eo 


to  satisfy 


eo  <  min 


(Qq  &i  b\  bi  \ 
c2  ’  4M,  ’  4Af,  ’  2c  J  ' 


* 

i 


Since  z(f)  »  F(z)(f)  it  follows  from  (18)  and  (23)  that  |z(<)|  <  2A/*«  whenever  c  <  <o. 
which  completes  the  proof. 
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6.  Open  Problem*.  There  ere  e  number  of  unresolved  mathematical  questions  which 
deserve  study  in  the  future.  These  are  the  following: 

(1)  Describe  the  dynamical  properties  of  flight  manuevers  in  the  vicinity  of  bifurcation 
points  on  the  equilibrium  set  E. 

(2)  Analyze  the  effect  of  small  random  disturbances  on  the  dynamics  of  flight  manuev¬ 
ers. 

(3)  Determine  what  effect  control  saturation  has  on  the  dynamics  of  aircraft  motion. 

(4)  Describe  the  dynamical  behavior  of  specific  flight  manuevers  with  inputs  H  with 
large  norm  N(H). 
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